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Abstract.  The  classical  maximum  flow  problem  sometimes  occurs  in  set¬ 
tings  in  which  the  arc  capacities  are  not  fixed  but  are  functions  of  a  single 
parameter,  and  the  goal  is  to  find  the  value  of  the  parameter  such  that 
the  corresponding  maximum  flow  or  minimum  cut  satisfies  some  side  con¬ 
dition.  Finding  the  desired  parameter  value  requires  solving  a  sequence 
of  related  maximum  flow  problems.  We  show  that  the  recent  maximum 
flow  algorithm  of  Goldberg  and  Tarjan  can  be  extended  to  solve  an  im¬ 
portant  class  of  such  parametric  maximum  flow  problems,  at  the  cost  of 
only  a  constant  factor  in  its  worst  case  time  bound.  Faster  algorithms  for 
a  variety  of  combinatorial  optimization  problems  follow  from  our  result. 

Keywords :  Algorithms,  data  structures,' graphs," maximum  flow,'  network 
flows;  networks. 

1.  Introduction 

The  well-known  maximum  flow  problem  calls  for  finding  a  maximum  flow  (or  alternatively 
a  minimum  cut)  in  a  capacitated  network.  This  problem  arises  in  a  variety  of  situations 
in  which  the  arc  capacities  are  not  fixed  but  are  functions  of  a  single  parameter,  and  the 
goal  is  to  find  the  value  of  the  parameter  such  that  the  corresponding  maximum  flow  (or 
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minimum  cut)  meets  some  side  condition.  The  usual  approach  to  solving  such  problems  is 
to  use  a  maximum  flow  algorithm  as  a  subroutine  and  use  either  binary  search,  monotonic 
search,  or  some  other  technique,  such  as  Megiddo’s  [28],  to  find  the  desired  value  of  the 
parameter. 

Existing  methods  take  no  advantage  of  the  similarity  of  the  successive  maximum  flow 
problems  that  must  be  solved.  In  this  paper,  we  address  the  question  of  whether  this 
similarity  can  lead  to  computational  efficiencies.  We  show  that  the  answer  to  this  question 
is  yes:  an  important  class  of  parametric  maximum  flow  problems  can  be  solved  by  extending 
the  new  maximum  flow  algorithm  devised  by  Goldberg  and  Tarjan  [13].  This  algorithm  is 
the  fastest  among  all  such  algorithms  for  real- valued  data,  uniformly  for  all  graph  densities. 
The  resulting  algorithm  for  the  parametric  problem  has  a  worst-case  time  bound  that  is 
only  a  constant  factor  greater  than  the  time  bound  to  solve  a  nonparametric  problem  of 
the  same  size.  The  parametric  problems  we  consider  are  those  in  which  the  capacities  of 
the  arcs  leaving  the  source  are  nondecreasing  functions  of  the  parameter,  those  of  arcs 
entering  the  sink  are  nonincreasing  functions  of  the  parameter,  and  those  of  all  other  arcs 
are  constant.  Our  parametric  maximum  flow  algorithm  has  a  variety  of  applications  in 
combinatorial  optimization. 

This  paper  consists  of  four  sections  in  addition  to  the  introduction.  In  Section  2  we  extend 
the  Goldberg-Tarjan  algorithm  to  find  maximum  flows  in  an  n- vertex,  m- arc  network  for 
O(n)  ordered  values  of  the  parameter  in  O(nmlog(n2/m))  time.  In  Section  3  we  use 
the  algorithm  of  Section  2  to  compute  information  about  the  minimum  cut  capacity  as 
a  function  of  the  parameter,  assuming  that  each  arc  capacity  is  a  linear  function  of  the 
parameter.  In  this  case,  the  minimum  cut  capacity  is  a  piecewise  linear  concave  function 
of  the  parameter.  We  describe  successively  more  complicated  algorithms  for  finding  the 
smallest  (or  largest)  breakpoint,  finding  a  maximum,  and  finding  all  the  breakpoints  of 
this  function.  Each  of  these  algorithms  runs  in  0(nmlog(nJ/m))  time. 

In  Section  4  we  discuss  applications  of  our  parametric  maximum  flow  algorithm  to  vari¬ 
ous  combinatorial  optimization  problems.  Depending  upon  the  application,  our  method 
is  faster  than  the  best  previously  known  method  by  a  factor  of  between  logn  and  n.  The 
applications  include  flow  sharing  problems  [3,  17,  20,  21,  26,  27],  zero-one  fractional  pro¬ 
gramming  problems  [4,  5,  11,  12,  23,  24,  29,  32,  33,  34,  38],  and  others  [9,  42].  Section  5 
contains  a  summary  of  our  results  and  some  final  remarks. 


2.  Parametric  maximum  flows  and  the  preflow  algorithm 


We  begin  in  this  section  by  reviewing  the  maximum  flow  algorithm  of  Goldberg  and  Tarjan 
[13],  here  called  the  preflow  algorithm.  Then  we  extend  their  method  to  the  parametric 
maximum  flow  problem  and  we  analyze  three  versions  of  the  parametric  preflow  algorithm. 
We  conclude  with  some  remarks  about  the  parametric  problem  and  our  algorithm  for 
solving  it. 

2.1.  Flow  terminology 

A  network  (see  [10,43])  is  a  directed  graph  G  =  (V,  F)  with  a  finite  vertex  set  V  and  arc  set 
E ,  having  a  distinguished  source  vertex  s ,  a  distinguished  sink  vertex  t,  and  a  nonnegative 
capacity  c(v,w)  for  each  arc  (v,w).  We  denote  the  number  of  vertices  by  n  and  the  number 
of  arcs  by  m.  We  assume  that  for  each  vertex  v ,  there  is  a  path  from  s  through  v  to  <;  this 
implies  n  =  O(m).  We  extend  the  capacity  function  to  arbitrary  vertex  pairs  by  defining 
c(v,w)  =  0  if  (u,u>)  E.  A  flow  f  on  G  is  a  nonnegative  function  on  vertex  pairs  satisfying 
the  following  three  constraints: 

f(v,w)  <  c(v,w)  for  (v,u>)  £  V  x  V  (capacity)  (2.1) 


f(v,w)  =  -f(w,v)  for  (v,tu)  E  V  x  V 
f(u,v)  =  0  for  v  £  V  —  {s,  t} 


(antisymmetry) 

(conservation) 


The  value  of  the  flow  /  is  5Z«i€V  /(v>0-  A  maximum  flow  is  a  flow  of  maximum  value. 

If  A  and  B  are  two  disjoint  vertex  subsets,  the  capacity  of  the  pair  A,B  is  c(A,B)  = 
52v£A,w£B  c(viw)-  A  cut  (X,X)  is  a  two-part  vertex  partition  (X  U  X  =  V,  X  n  X  =  0) 
such  that  s  £  X  and  i  £  X.  A  minimum  cut  is  a  cut  of  minimum  capacity.  If  /  is  a 
flow,  the  flow  across  the  cut  (X,X)  is  /( X,X)  =  f(v’w)-  The  conservation 

constraint  implies  that  the  flow  across  any  cut  is  equal  to  the  flow  value.  The  capacity 
constraint  implies  that  for  any  flow  /  and  any  cut  (AT,  X),  we  have  f(X,  X)  <  c(X,X), 
which  in  turn  implies  that  the  value  of  a  maximum  flow  is  no  greater  than  the  capacity 
of  a  minimum  cut.  The  max-flow  min-cut  theorem  of  Ford  and  Fulkerson  [10]  states  that 
these  two  quantities  are  equal. 


2.2.  The  preflow  algorithm 


The  preflow  algorithm  computes  a  maximum  flow  in  a  given  network.  To  describe  the 
algorithm  we  need  two  additional  concepts,  those  of  a  preflow  and  a  valid  labeling. 

A  preflow  f  on  G  is  a  real- valued  function  on  vertex  pairs  satisfying  the  capacity  constraint 
(2.1),  the  antisymmetry  constraint  (2.2),  and  the  following  relaxation  of  the  conservation 
constraint  (2.3): 

^  f(u,v )  >  0  for  all  v  £  V  -  {s}  (nonnegativity).  (2.4) 

u€V 

For  a  given  preflow,  we  define  the  excess  e(v)  of  a  vertex  v  to  be  /(u>v)  if  v  ^  s, 

or  infinity  if  v  =  s.  The  value  of  the  preflow  is  e(t).  We  call  a  vertex  r  ^  {s,<}  active  if 
e(r)  >  0.  A  preflow  is  a  flow  if  and  only  if  (2.4)  is  satisfied  with  equality  for  all  v  ^  {«,  t} , 
i.e.,  e(v)  =  0  for  all  v  ^  {s,  t}.  A  vertex  pair  (t >,w)  is  a  residual  arc  for  /  if  f{v,w )  <  c(v,u>); 
the  difference  c(v,w )  —  f(v,w)  is  the  residual  capacity  of  the  arc.  A  pair  (v,w)  that  is  not 
a  residual  arc  is  saturated.  A  path  of  residual  arcs  is  a  residual  path. 


A  valid  labeling  d  for  a  preflow  /  is  a  function  from  the  vertices  to  the  nonnegative  integers 
and  infinity,  such  that  d(t)  =  0,  d(s)  =  n,  and  d(v)  <  d( w)  +  l  for  every  residual  arc  (v,tu). 
The  residual  distance  dj(v,  w)  from  a  vertex  v  to  a  vertex  w  is  the  minimum  number  of  arcs 
on  a  residual  path  from  v  to  w,  or  infinity  if  there  is  no  such  path.  A  proof  by  induction 
shows  that  if  d  is  a  valid  labeling,  d(v)  <  min {df(v,t),df(v,s)  -f  n}  for  any  vertex  v. 


The  preflow  algorithm  maintains  a  preflow  /,  initially  equal  to  the  arc  capacities  on  arcs 
leaving  s  and  zero  on  other  arcs.  It  improves  /  by  pushing  flow  excess  toward  the  sink 
along  arcs  estimated  (by  using  d)  to  be  on  shortest  residual  paths.  The  value  of  /  gradually 
becomes  larger,  and  /  eventually  becomes  a  flow  of  maximum  value.  As  a  distance  estimate, 
the  algorithm  uses  a  valid  labeling  d ,  initially  defined  by  d(s)  =  n,d(v )  =  0  for  v  /  s.  This 
labeling  increases  as  flow  excess  is  moved  among  vertices;  such  movement  causes  residual 
arcs  to  change. 

To  implement  this  approach,  the  algorithm  uses  an  incidence  list  I(v)  for  each  vertex  v. 
The  elements  of  such  a  list,  called  edges,  are  the  unordered  pairs  {r,  u>}  such  that  (r,  w)  £  E 
or  (tL-,v)  £  E.  Of  the  edges  on  I(v),  one,  initially  the  first,  is  designated  the  current  edge 
of  v.  The  incidence  lists  I(v)  for  all  v  £  V  can  be  generated  from  an  arbirarily-ordered  arc 
list  E  in  O(m)  time. 


The  algorithm  consists  of  repeating  the  following  step  until  there  are  no  active  vertices: 

Push/Relabel.  Select  any  active  vertex  v.  Let  {v,u>}  be  the  current  edge  of  v.  Apply 
the  appropriate  one  of  the  following  three  cases: 

Push.  If  d(v)  >  d(w )  and  f(v,w )  <  c(v,w ),  send  6  =  min{e(v),c(r,tn)  -  /(v,tn)} 
units  of  flow  from  v  to  w.  This  is  done  by  increasing  f(v,w)  and  e(u>)  by  6 ,  and  by 
decreasing  f(w,v)  and  e(v)  by  6.  (The  push  is  saturating  if  6  =  c(v,w)  —  f(v,w)  and 
nonsaturating  otherwise.) 

Get  Next  Edge.  If  d(v)  <  d(w)  or  f(v,w)  =  c(v,u>),  and  {v,tu}  is  not  the  last  edge  on 
J(r),  replace  {v,tn}  as  the  current  edge  of  v  by  the  next  edge  on  I(v). 

Relabel  If  d(v)  <  d(w)  or  f(v,  w)  —  c(v,u>),  and  is  the  last  edge  on  J(t>), 

replace  d(v)  by  min{<f(tn)  |  {v,tn}  €  J(v)  and  f(v,w)  <  c(v,u>)}  +  1  and  make  the  first 
edge  on  I(v)  the  current  edge  of  v. 

When  the  algorithm  terminates,  /  is  a  maximum  flow.  A  minimum  cut  can  be  computed 
as  follows.  For  each  vertex  v,  replace  d(v )  by  min{rfj(v,s)  -f  n,  df(v,t )}  for  each  v  £  V. 
(This  replacement  cannot  decrease  any  distance  label.)  Then,  the  cut  (A',  A')  defined  by 
A'  =  {r  |d(v)  >  n}  is  a  minimum  cut  whose  sink  side  X  is  of  minimum  size,  a  property 
that  follows  from  Theorem  5.5  of  Ford  and  Fulkerson  [10]. 

The  efficiency  of  this  “generic”  form  of  the  preflow  algorithm  depends  upon  the  order  in 
which  active  vertices  are  selected  for  push/relabel  steps.  We  shall  consider  this  selection 
issue  after  we  have  extended  the  algorithm  to  the  parametric  problem.  For  the  moment, 
we  merely  note  the  bounds  derived  by  Goldberg  and  Tarjan  for  the  generic  algorithm  (with 
any  selection  order). 

Lemma  £.1  [13].  Any  active  vertex  v  has  dj(v,s)  <  oc,  which  implies  d(v)  <  2n  —  1.  The 
value  of  d(v)  never  decreases  during  the  running  of  the  algorithm.  The  total  number  of 
relabel  steps  is  thus  0(n2);  together  they  and  all  the  get  next  edge  steps  take  O(nm)  time. 

Lemma  £.£  [13].  The  number  of  saturating  push  steps  through  any  particular  residual  arc 
(t>,u>)  is  at  most  one  per  value  of  d(v).  The  total  number  of  saturating  push  steps  is  thus 
0(nm);  each  such  step  takes  0(1)  time. 


Lemma  t.S  [13],  The  total  number  of  nonsaturating  push  steps  is  0(n2m);  each  such  step 
takes  0(1)  time. 

In  all  variants  of  the  algorithm,  the  running  time  is  0(nm)  plus  0(1)  time  per  nonsaturat¬ 
ing  push  step;  making  the  algorithm  more  efficient  requires  reducing  the  number  of  such 
steps.  This  is  also  true  in  the  parametric  extension,  as  we  shall  see. 

2.3.  Extension  to  parametric  networks 

In  a  parametric  network ,  the  arc  capacities  are  functions  of  a  real-valued  parameter  A.  We 
denote  the  capacity  function  by  cA  and  make  the  following  assumptions: 

i.  cA(s,r)  is  a  nondecreasing  function  of  A  for  all  v  ^  t. 

ii.  cA(t’,<)  is  a  nonincreasing  function  of  A  for  all  v  /  s. 

iii.  cx(v,w)  is  constant  for  all  v  ^  s,  w  /  t. 

When  speaking  of  a  maximum  flow  or  minimum  cut  in  a  parametric  network,  we  mean 
maximum  or  minimum  for  some  particular  value  of  the  parameter  A. 

We  shall  address  the  problem  of  computing  maximum  flows  (or  minimum  cuts)  for  each 
member  of  an  increasing  sequence  of  parameter  values  Aj  <  A2  <  •  •  •  <  A;.  Successive 
values  are  given  on-line ;  that  is,  A|+1  need  not  be  known  until  after  the  maximum  flow 
for  A,  has  been  computed.  In  stating  time  bounds  we  shall  assume  that  the  capacity  of 
an  arc  can  be  computed  in  constant  time  given  the  value  of  A.  (Such  is  the  case  if,  for 
example,  the  arc  capacities  are  linear  functions  of  A.)  We  shall  also  assume  that  l  —  0(n); 
the  algorithm  we  shall  describe  computes  no  more  than  n  -  1  distinct  minimum  cuts,  no 
matter  how  many  values  of  A  are  given. 

We  shall  now  extend  the  preflow  algorithm  to  the  parametric  maximum  flow  problem. 
Suppose  that  for  some  value  A,  of  the  parameter  we  have  computed  a  maximum  flow  / 
and  a  valid  labeling  d  for  /.  What  is  the  effect  of  changing  the  value  of  the  parameter  to 
A^j?  The  capacity  of  each  arc  (s,v)  may  increase;  that  of  each  arc  (v,<)  may  decrease. 
Suppose  we  modify  /  by  replacing  f{v,t )  by  min{cA  (r,<),/(v,<)}  for  each  arc  (r,t)  £  E, 
and  replacing  /(s,v)  by  max{cA  (s,v),/(a,v)}  for  each  arc  (s,r)  £  E  such  that  d(v)  <  n. 
The  modified  /  is  a  preflow,  since  e(v)  for  v  {s,f}  can  only  have  increased.  Furthermore, 
d  is  a  valid  labeling  for  the  modified  /,  since  the  only  new  residual  arcs  are  of  the  form 


t 

I 

(j,t>)  for  d(v)  >  n  and  ( v,s )  for  d(v)  <  n.  This  means  that  we  can  compute  a  maximum 
flow  and  a  minimum  cut  for  AJ+1  by  applying  the  preflow  algorithm  beginning  with  the 
modified  /  and  the  current  d. 

This  idea  leads  to  the  following  parametric  preflow  algorithm.  The  algorithm  finds  a 

t 

maximum  flow  fx  and  a  minimum  cut  for  each  value  Xi  of  the  parameter.  It  ' 

consists  of  initializing  /  =  0,  d(s)  =  n,  d( v)  =  0  for  c^s,  and  i  =  0,  and  repeating  the  I 

following  three  steps  l  times:  , 

t 

Step  1.  (Update  preflow.)  Replace  i  by  i  +  1.  For  (v,t)  G  E,  replace  f{v,t)  by 
min{cA.(v,<),  f(v,t)}.  For  (s,t>)  G  E  with  d(v)  <  n,  replace  f(s,v)  by  max{cA.(s, v), 

/(*»«)}•  ; 

i 

Step  2.  (Find  maximum  flow.)  Apply  the  preflow  algorithm  to  the  network  with  arc 
capacities  corresponding  to  A{,  beginning  with  the  current  /  and  d.  Let  /  and  d  be 
the  resulting  flow  and  final  valid  labeling. 

Step  3.  (Find  minimum  cut.)  Redefine  d(y)  =  min{d^(r,s)  +  n,  dj(r,f)}  for  each 
v  €  V.  The  cut  (A\,  A'-)  is  then  given  by  A^  =  {v  \  <f(r)  >  n}. 

The  minimum  cuts  produced  by  the  algorithm  have  a  nesting  property  that  was  previously 
observed  in  the  context  of  various  applications  by  Eisner  and  Severance  [9],  Stone  [42], 
and  perhaps  others.  Megiddo  [26]  has  also  noted  a  similar  property  in  a  related  problem. 

Here,  the  result  follows  directly  from  our  algorithm. 

Lemma  2.4-  For  a  given  on-line  sequence  of  parameter  values  A}  <  A2  <  •••  <  AJf 
the  parametric  preflow  algorithm  correctly  computes  maximum  flows  ,  /2 ,  -  -  - ,  and 
minimum  cuts  (A'5 ,  A'j ), (A"2. X2),  •  •  • , (A”|,  Arj)  such  that  A'j  CX2  C  ■••C  A'(. 

Proof.  The  correctness  of  the  algorithm  is  immediate.  For  any  vertex  r,  the  label  d(v) 
cannot  decrease  in  Step  3,  which  implies  that  d(v)  never  decreases  during  the  running  of 
the  algorithm.  This  means  that  X}  C  X2  C  •  •  •  C  Xt,  which  in  turn  implies  that  there  can 
be  at  most  n  -  1  distinct  sets  Xx.  □ 
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2.4.  Analysis  of  the  parametric  preflow  algorithm 

In  view  of  Lemma  2.4,  Lemmas  2.1  and  2.2  hold  without  change  for  the  parametric  preflow 
algorithm.  Furthermore,  the  time  spent  in  Steps  1  and  3  is  0(m)  per  iteration  for  a  total 
of  0(lm )  =  O(nm)  time  (assuming  l  =  0(n)).  Thus  the  parametric  preflow  algorithm 
runs  in  O(nm)  time  plus  0(1)  time  per  nonsaturating  push. 

The  number  of  nonsaturating  pushes  depends  upon  the  order  in  which  push/relabel  steps 
are  performed.  We  shall  analyze  three  versions  of  the  parametric  preflow  algorithm.  For 
each,  we  show  that  the  time  bound  for  the  nonparametric  case  extends  to  the  parametric 
case  with  an  increase  of  at  most  a  constant  factor.  The  proofs  of  the  following  three 
theorems  are  analogous  to  those  given  in  [13]. 

We  first  analyze  the  generic  version,  in  which  push/relabel  steps  take  place  in  any  order. 

Theorem  2.5.  For  any  order  of  push/relabel  step  selection,  the  number  of  nonsaturating 
push  steps,  and  hence  the  running  time  of  the  parametric  preflow  algorithm,  is  0(n2m). 

Proof.  Let  =  £]{d(r)  I  v  is  active}  if  some  vertex  is  active,  and  $  =  0  otherwise.  A 
nonsaturating  push  step  decreases  $  by  at  least  one.  The  function  $  is  always  in  the  range 
0  to  2n2.  Step  1  can  increase  $  by  at  most  2n2,  for  a  total  over  all  iterations  of  Step  1  of 
0(ln2)  =  0(n3).  A  relabeling  step  increases  $  by  the  amount  by  which  the  label  changes. 
Thus  the  total  increase  in  $  due  to  relabeling  steps  is  0(n2).  A  saturating  push  step  can 
increase  $  by  at  most  2n.  Thus  the  total  increase  in  $  due  to  such  steps  is  0(n2m).  These 
are  the  only  ways  in  which  $  can  increase.  The  total  number  of  nonsaturating  push  steps 
is  bounded  by  the  total  increase  in  $  over  the  algorithm,  which  is  0(n2m).  □ 

Next  we  consider  the  first-in  first-out  (FIFO)  version  of  the  preflow  algorithm,  which  solves 
the  nonparametric  problem  in  0(n3)  time  [13:.  In  this  version,  a  queue  Q  is  used  to  select 
vertices  for  push/relabel  steps.  Initially  Q  is  empty.  At  the  beginning  of  Step  2  of  the 
parametric  preflow  algorithm,  every  active  vertex  is  appended  to  Q.  The  FIFO  algorithm 
consists  of  repeating  the  following  step  until  Q  is  empty: 

Discharge.  Remove  the  vertex  v  on  the  front  of  Q.  Apply  push/relabel  steps  to  r  until 
v  is  no  longer  active  or  v  is  relabeled.  If  a  push  from  v  to  another  vertex  w  makes  iv 
active,  add  w  to  the  rear  of  Q. 


Theorem  2.6.  For  the  FIFO  version  of  the  parametric  preflow  algorithm,  the  number  of 
nonsaturating  push  steps,  and  hence  the  total  running  time,  is  0(ns). 

Proof.  We  define  passes  over  the  queue  Q  as  follows.  The  first  pass  during  an  iteration 
of  Step  2  consists  of  the  discharge  steps  applied  to  the  vertices  initially  on  Q.  Each  pass 
after  the  first  in  an  iteration  consists  of  the  discharge  steps  applied  to  the  vertices  added 
to  Q  during  the  previous  pass.  There  is  at  most  one  nonsaturating  push  step  per  vertex  v 
per  pass,  since  such  a  step  reduces  e(v)  to  zero.  We  claim  that  the  total  number  of  passes 
is  0(n2  -f  In),  from  •which  the  theorem  follows. 

To  establish  the  claim,  we  define  #  =  max{d(t;)  |v  €  (?)  if  Q  ^  0,  and  $  =  0  if  Q  =  0. 
Consider  the  effect  on  $  of  a  pass  over  Q.  If  v  £  Q  at  the  beginning  of  a  pass  and  d(v)  =  $, 
then  v  $  Q  at  the  end  of  the  pass  unless  a  relabel  step  occurs  during  the  pass.  Thus,  if 
#  is  the  same  at  the  end  as  at  the  beginning  of  the  pass,  some  vertex  label  must  have 
increased  by  at  least  one.  If  $  increases  over  the  pass,  some  d(v)  must  increase  by  at  least 
the  amount  of  the  increase  in  4>.  From  the  end  of  one  iteration  to  the  beginning  of  the 
next,  can  increase  by  O(n).  Thus  (i)  the  total  number  of  passes  in  which  can  increase 
or  stay  the  same  is  0(n2  -f  In);  (ii)  the  total  number  of  passes  in  which  4>  can  decrease  is  at 
most  the  total  increase  in  $  between  passes  and  during  passes  in  which  it  increases,  which 
is  also  0(n2  -f  In).  We  conclude  that  the  total  number  of  passes  is  at  most  0(n2  4-  In), 
verifying  the  claim  and  hence  the  theorem.  □ 

A  more  elaborate  version  of  the  preflow  algorithm  [13]  uses  the  dynamic  tree  data  structure 
of  Sleator  and  Tarjan  [39,40  to  reduce  the  running  time  to  0(nm  log(n2/m))  .  The  cor¬ 
responding  version  of  the  parametric  flow  algorithm  also  runs  in  0(nm  log(n2/m))  time. 
It  uses  a  queue  Q  for  vertex  selection,  and  it  performs  discharge  steps  exactly  as  does 
the  FIFO  algorithm,  but  in  place  of  push/relabel  steps  it  uses  more  complicated  tree- 
push/relabel  steps.  A  tree-push/relabel  step  can  move  an  amount  of  flow  excess  through 
several  arcs  at  once.  Extending  the  analysis  in  [13]  to  the  parametric  case  is  straightfor¬ 
ward.  We  shall  merely  summarize  the  results. 

The  dynamic  tree  algorithm  uses  a  parameter  k,  the  maximum  tree  size,  which  can  be 
chosen  freely  in  the  range  from  2  to  n.  An  easy  extension  of  the  analysis  in  13  shows 
that  the  parametric  preflow  algorithm  runs  in  0(nm\ogk)  time  plus  O(logk)  time  per 
addition  of  a  vertex  to  Q.  Furthermore,  the  number  of  additions  of  a  vertex  to  Q  is  O(nm) 


plus  0(n/k)  per  pass  over  Q,  where  passes  are  defined  as  in  the  proof  of  Theorem  2.6. 
The  0(n2)  bound  on  the  number  of  passes  in  the  proof  of  Theorem  2.6  remains  valid 
if  the  dynamic  tree  algorithm  is  used  in  place  of  the  FIFO  algorithm.  Combining  these 
estimates,  we  obtain  an  0((nm  +  n3 /fc)logfc)  bound  on  the  total  running  time.  Choosing 
Jfc  =  max{2,  n2/m}  gives  the  following  theorem: 

Theorem  2.7.  The  parametric  preflow  algorithm  implemented  using  dynamic  trees  runs  in 
0(nmlog(n2/m))  time. 

2.5.  Additional  observations 

We  conclude  this  section  with  several  observations  about  the  parametric  maximum  flow 
problem  and  our  algorithm  for  solving  it.  Our  first  observation  concerns  variants  of  the 
parametric  maximum  flow  problem.  Our  algorithm  remains  valid  if  the  arc  capacity  func¬ 
tions  are  nonincreasing  on  arcs  out  of  s  and  nondecreasing  on  arcs  into  t,  and  the  values 
of  the  parameter  A  are  given  in  decreasing  order.  To  see  this,  merely  substitute  —  A  for  A. 
The  algorithm  also  applies  if  the  arc  capacity  functions  are  nondecreasing  on  arcs  out  of  s 
and  nonincreasing  on  arcs  into  i,  and  the  values  of  A  are  given  in  decreasing  order.  In  this 
case,  reverse  the  directions  of  all  the  arcs,  exchange  the  source  and  sink,  and  apply  the 
original  algorithm  to  this  reversed  network,  which  we  shall  denote  by  GR.  Each  minimum 
cut  (A",  A”)  generated  for  GR  will  correspond  to  a  minimum  cut  (A',  A”)  in  the  original 
(nonreversed)  network  that  will  have  the  source  side  of  minimum  size  instead  of  the  sink 
side.  Successively  generated  minimum  cuts  in  GR  will  correspond  to  cuts  with  successively 
smaller  source  sides  in  the  original  network. 

If  we  are  only  interested  in  computing  minimum  cuts,  there  is  a  variant  of  the  preflow 
algorithm  that  does  less  computation,  although  it  has  the  same  asymptotic  time  bound  [13. 
This  variant,  here  called  the  min-cut  preflow  algorithm ,  computes  a  preflow  of  maximum 
value  and  a  minimum  cut,  but  not  a  maximum  flow.  It  differs  from  the  original  algorithm 
in  that  a  vertex  v  is  considered  to  be  active  only  if  e(r)  >  0  and  d{ v)  <  n.  The  algorithm 
terminates  having  computed  a  maximum  preflow.  (A  preflow  /  is  maximum  if  and  only 
if  for  every  vertex  v,  v  ^  t  or  e(u)  =  0  or  dj(v,t)  <  oc.)  If  this  variant  is  used  to 
compute  minimum  cuts,  a  maximum  flow  for  a  desired  parameter  value  can  be  computed  by 
beginning  with  the  corresponding  maximum  preflow  and  converting  it  into  a  maximum  flow 
using  the  original  preflow  algorithm.  Most  of  the  applications  we  shall  consider  only  require 
the  computation  of  minimum  cuts  or  even  minimum  cut  values  and  not  maximum  flows. 
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So  far  we  have  required  all  arc  capacities  to  be  nonnegative,  but  if  we  are  only  interested 
in  computing  minimum  cuts,  we  can  allow  negative  capacities  on  arcs  out  of  a  and  on  arcs 
into  t.  This  is  because  there  is  a  simple  transformation  that  makes  such  arc  capacities 
nonnegative  without  affecting  minimum  cuts  {34].  For  a  given  vertex  v,  suppose  we  add 
a  constant  A(t>)  to  c(a,v)  and  c(t>,<).  Then  the  minimum  cuts  do  not  change  since  the 
capacity  of  every  cut  is  increased  by  A(u). 

By  adding  a  suitably  large  A(r)  to  c(a,u)  and  c(v,t)  for  each  v,  we  can  make  all  the  arc 
capacities  positive.  In  the  parametric  problem,  we  can  choose  a  new  function  A  on  the 
vertices  of  G  for  each  new  value  A^  of  A  without  affecting  the  0(nmlog(n2/m))  time  bound 
of  our  algorithm.  It  suffices  to  choose 

AAj(v)  =  max{0,cAi(a,r)}  +  max{0, -cA](t>,<)}, 

=  aa»  +  cx.(v,t)  -  cA.+l(r,<)  for  i  >  1. 

With  this  choice,  the  transformed  arc  capacities  are  nondecreasing  functions  of  A  on  arcs 
leaving  a  and  constant  on  arcs  that  enter  t.  Although  the  same  effect  could  be  obtained 
by  adding  a  sufficiently  large  constant  to  the  capacities  of  these  arcs,  the  modification 
we  have  described  has  the  additional  advantage  of  keeping  capacities  as  small  as  possible. 
In  subsequent  sections,  when  discussing  minimum  cut  problems,  we  shall  allow  arbitrary 
capacities,  positive  or  negative,  on  arcs  leaving  a  and  arcs  entering  t. 

The  minimum  cuts  corresponding  to  various  parameter  values  have  a  nesting  property  that 
is  a  strengthening  of  Lemma  2.4.  The  following  lemma  is  an  extension  of  known  results 
[10,  p.13]  that  we  shall  need  in  the  next  section.  To  prove  the  lemma,  we  shall  use  the 
min-cut  preflow  algorithm  discussed  above. 

Lemma  2.8.  Let  ( X,X )  be  any  minimum  cut  for  A  =  Aj,  and  let  (F,  F)  be  any  minimum 
cut  for  A  =  A2  such  that  Aj  <  A2.  Then,  (AT  D  F, X  U  F)  is  a  minimum  cut  for  A  =  Aj  and 
(X  U  F,  X  fl  F)  is  a  minimum  cut  for  A  =  A2. 

Proof  Run  the  min-cut  parametric  algorithm  for  A  =  A 2 ,  followed  by  A  =  A2.  At  the 
beginning  of  the  computation  for  A  =  A2,  all  vertices  v  €  X  -  {a}  have  d(v)  >  n.  Thus, 
after  a  maximum  preflow  is  computed  for  A  =  A2,  all  arcs  (v,w)  with  v  €  X,w  £  A'  are 
saturated  (their  flow  has  not  changed  during  the  computation  for  A  =  A2.)  Since  (}', F)  is 
a  minimum  cut  for  A  =  A2,  all  arcs  (v,w)  with  v  €  F,  w  €  F  are  saturated.  Furthermore, 


if  v  6  Y  -  {*}  then  e(v)  =  0,  since  the  net  flow  across  (Y,  F)  must  be  equal  to  the  excess 
at  t. 

Now  consider  the  cut  Z  =  (Ar  U  Y ),  2  =  (X  n  F).  Any  arc  (v,w)  w’ith  v  £  Z ,  w  £  Z  must 
be  staurated.  Since  v  £  Z  —  {t}  implies  e(v)  =  0,  the  cut  ( Z,Z )  must  have  capacity  e(t), 
and  hence  it  must  be  a  minimum  cut. 

The  proof  for  (A”  D  Y,  X  U  F)  is  similar:  proceed  on  GR,  and  run  the  min-cut  parametric 
algorithm  for  A  =  A2,  followed  by  A  =  Aj.  □ 

A  direct  consequence  of  this  lemma  is  the  following  corollary,  which  we  shall  need  in  the 
next  section. 

Corollary  2.9.  Let  (Arj,Aj)  be  a  minimum  cut  for  A  =  A1?  let  (Ar2,A'2)  be  a  minimum  cut 
for  A  =  A2  such  that  Xj  C  X2  and  A3  <  A2,  and  let  A3  be  such  that  Aj  <  A3  <  A2.  Then 
there  is  a  cut  (X3 ,  X3 )  minimum  for  A3  such  that  Arj  C  A'3  C  X2. 

Proof.  Let  (Ar3,A'3)  be  any  minimum  cut  for  A  =  A3.  Take  Ar3  =  (X'3  U  A'j)  n  A'2, 
A'3  =  V’  -  X3,  and  apply  Lemma  2.8  twice.  □ 

Our  last  observation  is  that  if  the  graph  G  is  bipartite,  the  0(nm  log(n2/m))  time  bound 
for  computing  parametric  maximum  flows  can  be  improved  slightly.  Suppose  V  =  A  U  B, 
A  n  B  ^  0,  and  every  arc  in  G  has  one  vertex  in  A  and  one  in  B.  Let  nA  =  |A|  and 
nB  =  |J3j  and  suppose  that  nA  <  nB.  Then  the  time  to  compute  maximum  flows  for  / 
ordered  values  of  A  can  be  reduced  to  0(nAm  log(nA/m  -f  2))  if  /  =  0(nA).  This  requires 
modifying  the  preflow  algorithm  so  that  only  vertices  in  A  are  active,  and  modifying  the 
use  of  the  dynamic  tree  data  structure  so  that  such  a  tree  contains  as  many  vertices  in  A  as 
in  B.  The  details  can  be  found  in  [41].  The  bound  is  slightly  worse  if  /  =  (see  [41]). 

3.  The  min-cut  capacity  function  of  a  parametric  network 

For  a  parametric  network,  we  define  the  min-cut  capacity  function  *(A)  to  be  the  capacity 
of  a  minimum  cut  as  a  function  of  the  parameter  A.  We  shall  assume  throughout  this 
section  that  the  arc  capacities  are  linear  functions  of  A  satisfying  the  conditions  (i)-(iii)  of 
Section  2.  It  is  well  known  [9,42]  and  follows  from  the  results  of  Section  2  that  under  this 
assumption  k(A)  is  a  piecewise-linear  concave  function  with  at  most  n  — 2  breakpoints.  (By 
a  breakpoint  we  mean  a  value  of  A  at  which  the  slope  of  *(A)  changes.)  The  n  -  1  or  fewer 
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line  segments  forming  the  graph  of  /c(A)  correspond  to  n  -  1  or  fewer  distinct  cuts.  We 
shall  develop  three  algorithms  for  computing  information  about  k(X).  The  first  computes 
the  smallest  (or  equivalently  the  largest)  breakpoint.  The  second  computes  a  value  of 
A  at  which  /c(A)  is  maximum.  The  third  computes  all  the  breakpoints.  Each  of  these 
algorithms  uses  the  algorithm  of  Section  2  as  a  subroutine  and  runs  in  0(nmlog(n2/m)) 
time.  Although  the  algorithm  for  computing  all  breakpoints  solves  all  three  problems,  we 
shall  present  all  three  algorithms  since  each  is  more  complicated  than  the  preceding  one 
and  since  the  resulting  difference  in  constant  factors  may  be  important  in  practice. 

We  shall  assume  that  the  capacities  cA(s, v)  and  cx(v,t)  are  given  in  the  form  cA(s,t>)  = 
a0(r)  +  Aaj(n)  and  cx(v,i)  =  b0(v )  -  Xb1(v),  with  arbitrary  coefficients  a0,60  and  nonneg¬ 
ative  coefficients  a1,b1.  A  minimum  cut  (A0,A0)  for  some  A  =  A0  gives  an  equation  for  a 
line  that  contributes  a  line  segment  to  the  function  k(A)  at  A  =  A0.  This  line  is  Lx o(A)  = 
o0  +  A 0Ot  where  a0  =  cAo(A0,A'0)  -  XO0O  and  0O  =  '£v£Xo  o,(i>)  -  Mv)- 

3.1.  Computing  the  smallest  breakpoint  of  /c(A) 

To  compute  the  smallest  breakpoint  of  /c(A)  we  use  an  algorithm  stated  by  Gusfield  [1 7]  for 
an  application  involving  scheduling  transmissions  in  a  communication  network,  discussed 


in  more  detail  in  Section  4.1. 


The  algorithm  consists  of  the  following  two  steps. 


Step  0.  Compute  Aj,A2  such  that  the  smallest  breakpoint  A0  satisfies  A2  <  A0  <  A2. 
Compute  a  cut  (A'j,  A, )  that  is  a  minimum  for  Ar  Go  to  Step  1. 


Step  1.  Compute  a  cut  (A2,  A"2)  that  is  a  minimum  for  A2.  If  LXi( Aj)  =  LXi(X2), 
stop:  A2  is  the  smallest  breakpoint.  Otherwise,  replace  A2  by  the  value  of  A  such 
that  LXx( A)  =  LXi(X)  and  repeat  Step  1.  (The  appropriate  value  of  A  is  (a2  -  a2)/ 

(ft  -  A)-) 


The  values  of  A2  generated  by  this  algorithm  are  strictly  decreasing;  thus  the  para¬ 
metric  preflow  algorithm  of  Section  2  applied  to  GR  performs  all  ierations  of  Step  1  in 
0(nmlog(n2/m))  time,  saving  a  factor  of  n  over  Gusfield’s  algorithm  [17]. 
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In  Step  0,  it  suffices  to  select  A3  sufficiently  small  that  for  each  vertex  v  such  that  (a,v) 
or  (v,t)  is  of  nonconstant  capacity,  cAj(s,t>)  4-  53u€y-{*  tj  e(u,t>)  <  cAj(v,<).  A  suitable 
value  of  Aj  is 


Ai  =  ,5“?  A  6o(v)  -  ao(v)  - 
t>€  V  — {#,t} 


/wEV-{i,l} 


c(u,v) 


aj(v)  +  6j(r)  >  0  }  -1.  (3.1) 


aj(u)  +  6j(u) 

Similarly,  it  suffices  to  select  A2  sufficiently  large  that  for  each  vertex  v  such  that  (s,v) 
or  (v,f)  is  of  nonconstant  capacity,  cXj(v,t)  +  c(v,w)  <  cXj(s,v).  A  suitable 


value  of  A2  is 


x  ft/,  ,  ,  ,  ‘K™) 

A,  =  max  {  Mv)  _  °o(v)  -f - r"T — r~7~^ 

2  v€V-{,,«}1  ov  '  ov  ’  Oj(«)  + 


a3(u)  +  ij(u)  >  0  }  +1.  (3.2) 


Essentially  the  same  algorithm  can  be  used  for  computing  the  largest  breakpoint;  instead 
of  successively  decreasing  A2  and  using  GR ,  successively  increase  A3  and  use  G. 


3.2.  Finding  a  maximum  of  /c(A) 


Our  algorithm  for  finding  the  value  of  A  that  maximizes  /c(A)  is  based  on  a  simple  method 
of  iterative  interval  contraction  for  computing  the  maximum  /(A*)  of  a  strictly  concave 
and  continuously  differentiable  function  /(A)  on  a  nonempty  interval  [A3,A2]  of  the  real 
line.  The  method  is  as  follows.  First,  compute  the  function  values  and  the  tangents  of  /(A) 
at  each  end  of  the  given  interval.  Second,  compute  the  point  Aj  €  [A a ,  A2 ]  where  the  two 
tangent  lines  intersect,  and  also  compute  /'( A3).  Then,  if  /'( A3)  <  0  replace  A2  by  A3  and 
repeat;  if  /'( A3)  >  0  replace  A3  by  A3  and  repeat;  if  /'( A3)  =  0  accept  A3  as  the  solution. 
Of  course  this  algorithm  need  not  terminate,  but  it  will  converge  to  the  maximum. 


The  method  is  seldom  used  in  this  general  setting  because  it  is  inferior  to  several  other 
algorithms  for  one- dimensional  maximization.  But  it  can  be  specialized  in  the  obvious  way 
to  handle  the  piecewise-linear  concave  function  k(A)  efficiently.  A  maximum  of  *c(A)  can 
be  computed  in  as  many  function  evaluations  as  there  are  linear  segments  that  comprise 
«( A),  namely  n—  1  or  fewer.  Using  the  notation  introduced  above,  \3  =  (q2 -q3)/(/?j  —02) 
if  the  line  segments  of  k(A)  at  As  and  at  A2  are  distinct.  Otherwise,  the  search  terminates 
with  a  line  segment  of  zero  slope  and  A*  =  A3  (or  A*  =  A2).  The  algorithm  will  compute  a 
maximum  of  *(A)  in  0(n2mlog(n2/m))  time  since  at  most  n  —  1  minimum  »ut  problems 
must  be  solved. 


The  running  time  of  this  algorithm  can  be  improved  by  partitioning  the  sequence  of  succes¬ 
sive  values  of  A3  into  two  subsequences,  one  increasing  and  the  other  decreasing.  It  is  then 


possible  to  use  two  concurrent  invocations  of  the  parametric  preflow  algorithm:  Invocation 
I  that  starts  with  A,  and  computes  minimum  cuts  of  G  for  an  increasing  sequence  of  A 
values,  and  Invocation  D  that  starts  tvith  A2  and  computes  minimum  cuts  of  GR  for  a 
decreasing  sequence  of  A  values.  A  new  value  A3  is  a  member  of  the  increasing  sequence  if 
>  0,  and  a  member  of  the  decreasing  sequence  otherwise.  The  initial  values  of  A,  and 
A2  must  be  such  that  all  breakpoints  lie  in  the  interval  [A3 ,  A2],  a  property  that  is  assured 
by  the  values  of  A  defined  by  (3.1)  and  (3.2). 

The  algorithm  to  compute  a  maximum  of  «( A)  consists  of  the  following  four  steps. 

Step  0.  Compute  the  initial  values  Aj  and  A2  from  (3.1)  and  (3.2).  Start  concurrent 
invocations  (I  and  D )  of  the  parametric  preflow  algorithm  of  Section  2:  For  A  =  An 
Invocation  I  computes  a  minimum  cut  (A',  ,Xt)  having  |_Vj  |  maximum;  for  A  =  A2, 
Invocation  D  computes  a  minimum  cut  (A2,  A'2)  having  |Ar2|  minimum. 

Step  1.  Compute  A3  =  (a2  —  o  ])/(/?]  -  02),  pass  A3  to  both  invocations  I  and 
D,  and  run  them  concurrently.  If  invocation  I  finds  a  minimum  cut  (A/,  A'/)  first, 
suspend  invocation  D  and  go  to  Step  2  (the  other  case  is  symmetric.)  Compute 

^3  =  ai(v)  ~  Mv)- 

Step  2.  If  03  =  0,  stop:  A*  =  A3.  Otherwise,  if  03  >  0,  replace  Aj  by  A3,  back  up 
invocation  D  to  its  state  before  it  began  processing  A3,  and  go  to  Step  1.  Otherwise, 
go  to  Step  3. 

Step  S  (03  <  0).  Finish  running  the  invocation  D  on  A3.  This  produces  a  minimum 
cut  (XfjJCf),  not  necessarily  the  same  as  ( X3,X3 ).  If  03  >  0,  stop:  A*  =  A3. 
Otherwise,  replace  A2  by  A3,  back  up  invocation  I  to  its  state  before  processing  A3. 
and  go  to  Step  1. 

Backing  up  invocation  U  or  /  as  required  in  Steps  2  and  3  is  merely  a  matter  of  restoring 
the  appropriate  flow  and  valid  labeling,  which  takes  0(m)  time.  The  total  time  spent 
during  one  iteration  of  Steps  1,  2  and  3  is  proportional  to  the  time  spent  in  invocation  I 
or  D,  whichever  one  is  run  to  completion  on  A3  and  not  backed  up.  The  total  number  of 
values  of  A3  processed  is  0(n).  Thus  the  total  time  is  proportional  to  the  time  necessary  to 
run  the  parametric  preflow  algorithm  of  Section  2.3  twice,  once  on  an  increasing  sequence 
of  values  and  once  on  a  decreasing  sequence  of  values;  i.e.,  O(nmlog(n2/m))  . 
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3.3.  Finding  all  breakpoints  of  «(A) 

In  some  applications  it  is  necessary  to  produce  all  the  line  segments  or  breakpoints  of  «(A), 
possibly  along  with  the  corresponding  minimum  cuts.  To  do  this  we  extend  the  maximum¬ 
finding  algorithm  of  the  previous  section.  This  algorithm  uses  iterative  contraction  of  the 
interval  [Aj,A2];  it  ignores  breakpoints  that  lie  in  the  discarded  portion  of  the  interval. 
We  can  find  all  the  breakpoints  by  proceeding  as  in  the  algorithm  of  Section  3.2  but 
using  a  divide-and-conquer  strategy  that  recursively  examines  both  of  the  subintervals 
[Aj,Aj]  and  [A3,A2]  into  which  the  current  interval  is  split  by  the  new  value  As.  This 
method  was  proposed  by  Eisner  and  Severance  [9]  for  bipartite  graphs  in  the  context  of  a 
database  record  segmentation  problem  (see  Section  4.4).  Unfortunately,  a  straighforward 
implementation  of  this  idea  yields  an  0(n2mlog(n2/m))-time  algorithm.  To  obtain  a 
better  bound  it  is  necessary  to  use  two  concurrent  invocations  of  the  parametric  preflow 
algorithm,  and  also  use  graph  contraction  so  that  recursive  invocations  of  the  method 
compute  cuts  on  smaller  and  smaller  graphs. 

If  G  is  a  network  and  X  is  a  set  of  vertices  such  that  exactly  one  of  s  and  i  is  in  X, 
we  define  G(Ar),  the  contraction  of  G  by  X,  to  be  the  network  formed  by  shrinking  the 
vertices  in  X  to  a  single  vertex,  eliminating  loops,  and  combining  multiple  arcs  by  adding 
their  capacities.  The  algorithm  we  present  reports  only  the  breakpoints  of  k(A),  although 
it  computes  cuts  corresponding  to  the  line  segments  of  the  graph  of  a(A).  If  the  actual 
cuts  are  needed,  they  can  either  be  saved  as  the  computation  proceeds  or  computed  in  a 
postprocessing  step  using  one  application  of  the  method  in  Section  2.3. 

Our  algorithm  uses  a  recursive  procedure  called  slice.  With  each  network  G  to  which  slice 
is  applied,  we  associate  four  pieces  of  information:  Two  values  of  A,  denoted  by  A,  and  A2, 
and  two  flows  and  /2,  such  that  /j  is  a  maximum  flow  for  A5 ,  /2  is  a  maximum  flow  for 
A2,  the  cut  ({•s},V  —  {s})  is  the  unique  minimum  cut  for  Aj,  the  cut  (V  -  {<},{/})  is  the 
unique  minimum  cut  for  A2,  and  Aj  <  A2.  The  initial  values  for  Aj  and  A2  are  computed 
from  (3.1)  and  (3.2)  as  before.  The  breakpoint  algorithm  consists  of  the  following  three 


steps. 


Step  1.  Compute  A,  according  to  (3.1).  Compute  a  maximum  flow  /,  and  minimum 
cut  (A’pA,)  for  Aj  such  that  |Xj|  is  maximum  by  applying  the  preflow  algorithm  to 
G.  Let  G'  =  G{ Aj). 
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Step  t.  Compute  A2  according  to  (3.2).  Compute  a  maximum  flow  /2  and  minimum 
cut  (X2,X2)  for  A2  such  that  |Aj|  is  minimum  by  applying  the  preflow  algorithm  to 
(G')r.  Let  G"  =  <?'(A'2). 


Step  S.  If  G"  contains  at  least  three  vertices,  let  /"  and  f2  be  the  flows  in  G" 
corresponding  to  /3  and  f2  respectively;  perform  slice  (G",  Aj,  A2,/",/2  ),  where  slice 
is  defined  as  follows: 

Procedure  slice(G ,  Aj ,  A2 ,  fy ,  f2 ). 

Step  Si.  Let  A3  be  the  value  of  A  such  that  cAg({-*},V  -  {a})  =  ca»(^  - 
(This  value  will  satisfy  A3  <  A3  <  A2.) 

Step  S2.  Run  the  preflow  algorithm  for  the  value  A3  on  G  starting  with  the  preflow  /J 
formed  by  increasing  f}  on  arcs  (,s,t>)  to  saturate  them  and  decreasing  /3  on  arcs  (v,  t) 
to  meet  the  capacity  constraints.  As  an  initial  valid  labeling,  use  d(v)  =  min{d^,(v,t), 
dj,  ( v ,  s)  +  n}.  Concurrently,  run  the  preflow  algorithm  for  the  value  A3  on  GR  starting 
with  the  preflow  f2  formed  by  increasing  f2  on  arcs  ( v,t )  to  saturate  them  and  de¬ 
creasing  f2  on  arcs  (s,r)  to  meet  the  capacity  constraints.  As  an  initial  valid  labeling, 
use  d(v)  =  min{rfj.(s,u),  rfy,  (t,  v)  -f  n}.  Stop  when  one  of  the  concurrent  applications 
stops,  having  computed  a  maximum  flow  /3 .  Suppose  the  preflow  algorithm  applied 
to  G  stops  first.  (The  other  case  is  symmetric.)  Find  the  minimum  cuts  (A'3,X3) 
and  (Ar3,A”3)  for  A3  such  that  |A”3|  is  minimum  and  |A”3|  is  maximum.  If  |A”3|  >  n/2, 
complete  the  execution  of  the  preflow  algorithm  on  GR  and  let  f3  be  the  resulting 
maximum  flow. 

Step  SS.  If  cx(X3,X3)  ^  eA( X3,X3)  for  some  A,  report  A3  as  a  breakpoint. 

Step  S4.  If  X3  ^  {s},  perform  slice  (G(A3),  A3 ,  A3,/j  ,/3).  If  X'3  ^  {t},  perform 
slice  (G(X3 ),  A3 ,  A2 , /3 , f2 ). 

The  correctness  of  this  algorithm  follows  from  Corollary  2.9.  Note  that  the  minimum  cuts 
computed  in  Step  S2  correspond  to  minimum  cuts  for  A3  in  the  original  network,  with  the 
correspondence  obtained  by  expanding  the  contracted  vertex  sets.  Since  each  vertex  of  G 
is  in  at  most  one  of  the  two  subproblems  in  Step  S4,  there  are  O(n)  invocations  of  slice. 


3.4.  Analysis  of  the  breakpoint  algorithm 

Two  ideas  underly  the  efficiency  of  the  breakpoint  algorithm.  To  explain  them,  we  need  to 
develop  a  framework  for  the  analysis  of  the  algorithm.  We  shall  charge  to  an  invocation  of 
slice  the  time  spent  in  the  invocation,  not  including  the  time  spent  in  nested  invocations. 
The  time  charged  to  one  invocation  is  then  O(m)  plus  the  time  spent  running  the  preflow 
algorithm  in  Step  S2.  Summing  O(m)  over  all  O(n)  invocations  of  slice  gives  a  bound  of 
0(nm).  It  remains  to  estimate  the  time  spent  running  the  preflow  algorithm. 

The  first  idea  contributing  to  the  speed  of  the  algorithm  is  that  the  results  of  Section  2 
allow  us  to  bound  the  time  of  a  sequence  of  preflow  algorithm  applications,  not  just  a 
single  one,  by  0(nmlog(n2/m))  .  That  is,  if  we  charge  this  much  time  for  an  invocation 
of  slice,  we  can  regard  certain  of  the  nested  invocations  as  being  free.  The  second  idea 
is  that  running  the  preflow  algorithm  concurrently  on  G  and  on  GR  allows  us  to  regard 
the  larger  of  the  nested  invocations  in  Step  S4  as  being  free.  This  leads  to  a  recurrence 
bounding  the  running  time  whose  solution  is  0(nmlog(n2/m))  . 

Consider  an  invocation  of  slice  G(Aj ,  A2 , /2 , /2).  Let  G}  =  G(  X3)  as  computed  in  Step 
S4,  and  let  G2  =  G(A"3);  let  n^mj  and  n2,m2  be  the  numbers  of  vertices  and  arcs  in 
G j  and  G2 ,  respectively.  We  regard  this  invocation  of  slice  as  being  a  continuation  of 
the  algorithm  of  Section  2.3  applied  to  G,  with  Aj  the  most  recently  processed  value  of 
A  and  /j  the  resulting  maximum  flow.  Simultaneously,  we  regard  the  invocation  as  being 
a  continuation  of  the  algorithm  of  Section  2  applied  to  GR,  with  A2  the  most  recently 
processed  value  of  A  and  f2  the  resulting  flow. 

With  this  interpretation  we  can  regard  the  preflow  algorithm  applications  in  Step  S2  as 
being  free,  but  if  |Ar3|  <  n/2  we  must  account  for  new  applications  of  the  Section  2 
algorithm  on  G(Ar3)  and  Gfi(A3),  and  otherwise  (i.e.,  |A'3|  <  n/2)  we  must  account  for 
new  applications  of  the  Section  2  algorithm  on  G( Ar3)  and  Gfl(A3).  Thus  we  obtain  the 
following  bound  on  the  time  spent  in  invocations  of  the  preflow  algorithm.  If  G  has  n 
vertices  and  m  arcs,  the  time  spent  in  such  invocations  during  the  computation  of  k(A)  is 
at  most  T(n,m)+  0(nmlog(n2/m))  ,  where  T(n,m)  is  defined  recursively  as  follows: 


T(n,  m)  = 


0  if  n  <  3; 

max{T(ni,m1)  +  T(n2,m2 )  +  0(n,mj  log(n2,/m,))  : 
n},n2  >  3; Hj  +  n2  <  n  4-  2; 

<  n2;mj,m2  >  l;mj  +  m2  <  m  1}  if  n  >  3. 


Remark.  In  this  analysis,  the  sequence  of  preflow  algorithm  invocations  associated  with 
a  particular  application  of  the  algorithm  of  Section  2.3  is  on  a  sequence  of  successively 
smaller  graphs,  but  the  analysis  in  Section  2.4  remains  valid.  □ 


To  solve  the  recurrence  for  T(n,m),  we  begin  by  simplifying  it.  Observe  that  T(n,m)  = 
0(Tj(n,m)log(n2/m)),  where  Tj(n,m)  is  defined  as  follows: 


r3(n,m)  = 


0  if  n  <  3  ; 

max{Tj(nj ,mj)  +  Tj(n2,m2)  +  n,mj  : 
njin2  -  3;  rij  +  n2  <  n  4-  2; 

n3  <  n2;  rrij,  m2  >  1;  mj  +  m2  <  m  +  1}  if  n  >  3. 


By  making  the  change  of  variables  n'  =  n  —  2,  m'  —  m  —  1,  we  obtain  Tj(n,m)  = 
T2(n  -  2,m  -  1),  where  T2  is  defined  as  follows: 


{0  if  n  <  1; 

max{r2(nj,m3)  +  T2(n2,m2)  +  rijmj  -f  2m1  +  n 3  -r  2 
rij  ,n2  <  1;  n,  +  n2  <  n; 

n j  <  n2;  mj,m2  >  0;  m3  +  m2  <  m}  if  n  >  1. 


Since  the  recurrence  for  T2(n,m)  can  be  unwound  at  most  n  -  1  times  before  all  branches 
are  terminal,  the  additive  term  “2m,  -f  n3  +  2”  contributes  at  most  (n  —  l)(2m  +  n  -v  2) 
to  T2(n,m).  That  is,  T2(n,m)  =  T3(n,m)  +  O(nm),  where  T3  is  defined  as  follows: 


r3(n,m) 


'  0  if  n  <  1; 

max  {T3(nj,mj)  +  T3(n2,m2)  +  rijm, 
<  Tij ,  n2  <1;  Tij  4-  n2  <  n; 

rij  <  n2;  m1,m2  >  0; 

.  77i  3  +  m  2  <  m}  if  n  >  1. 


An  easy  proof  by  induction  shows  that  T3(7i,m)  <  nm.  This  implies  that  T(n,m)  is 
0(nrnlog(n2 /m))  and  in  turn  that  the  total  running  time  of  the  algorithm  is  of  the  same 
order. 


3.5.  Additional  observations 


We  conclude  this  section  with  two  observations.  First,  note  that  a  complete  set  of  minimum 
cuts  for  all  values  of  A  can  be  represented  in  O(n)  space:  store  with  each  vertex  v  &  {s,f} 
the  breakpoint  at  which  v  moves  from  the  sink  side  to  the  source  side  of  a  minimum  cut, 
for  a  set  of  minimum  cuts  whose  source  sides  are  nested.  The  breakpoint  algorithm  can 
be  augmented  to  compute  this  information  without  affecting  its  asymptotic  time  bound. 
Second,  the  time  bound  of  the  three  algorithms  in  Sections  3. 1-3.3  can  be  improved  to 
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0(nAm  log(n^/m  +  2))  if  G  is  bipartite  and  ic(A)  has  ^(n^)  breakpoints.  Here  nA  is 
the  size  of  the  smaller  half  of  the  bipartite  partition  of  V.  This  bound  follows  using  the 
bipartite  variant  of  the  preflow  algorithm  mentioned  at  the  end  of  Section  2.5  (see  [41]  for 
details). 

4.  Applications 

In  this  section,  we  give  a  number  of  applications  of  the  algorithms  in  Sections  2  and  3.  For 
each  application,  we  obtain  an  algorithm  running  in  0(nmlog(n2/m))  time,  where  n  is 
the  number  of  vertices  and  m  is  the  number  of  arcs  in  the  graph  involved  in  the  problem. 
For  applications  in  which  the  graph  is  bipartite,  the  bound  is  0(nAmlog(n2A/m  +  2)), 
where  nA  is  the  size  of  the  smaller  half  of  the  bipartite  partition  of  the  vertex  set.  (When 
the  latter  bound  is  applicable,  we  shall  state  it  within  square  brackets.)  Depending  upon 
the  application,  our  bound  is  a  factor  of  from  logn  to  n  better  than  the  best  previously 
known  bound.  Our  applications  fall  into  four  general  categories:  Flow  sharing  problems, 
zero-one  fractional  programming  problems,  parametric  zero-one  polynomial  programming 
problems,  and  miscellaneous  applications. 

4.1  Flow  sharing 

Consider  a  network  with  a  set  of  sources  S  =  {sj,s2,  ■  •  •  ,  st}  and  a  single  sink  t ,  in  which 
we  want  to  find  a  flow  from  the  sources  in  5  to  t.  We  require  flow  conservation  at  vertices 
not  in  5  o  {<}.  We  can  model  this  problem  as  an  ordinary  one-source,  one-sink  problem 
by  adding  a  supersource  s  and  an  arc  (s,s,)  of  infinite  capacity  for  each  i  £  {1,  •••,&}. 
The  resulting  network  can  have  many  different  maximum  flows,  with  different  utilizations 
of  the  various  sources;  we  define  the  utilization  ui  of  source  s,  to  be  the  flow  through  the 
arc  (5,5,).  The  question  arises  of  how  to  compare  the  quality  of  such  flows.  Suppose  each 
source  s,  has  a  positive  weight  u>,.  Severed  figures  of  merit  have  been  proposed,  leading  to 
the  following  optimization  problems: 

i.  Perfect  sharing :  Among  flows  with  ui/wx  equal  for  all  i  £  {1,  •  •  • ,  &},  find  one  that 
maximizes  the  flow  vedue  e(<). 

ii.  Mazimin  sharing:  Among  maximum  flows,  find  one  that  maximizes  the  smallest 

u./u-',,  *  €  -i*}- 


race 
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iii.  Minimax  sharing:  Among  maximum  flows,  find  one  that  minimizes  the  largest 

i  €  {l,--- 

iv.  Optimal  sharing:  Among  maximum  flows,  find  one  that  simultaneously  maximizes 
the  smallest  uijwi  and  minimizes  the  largest  ujw^  i  £  {1,  •••,&}. 

v.  Lexicographic  sharing:  Among  maximum  flows,  find  one  that  lexicographically 
maximizes  the  fc-component  vector  whose  j-th  component  is  the  y-th  smallest  ui/wi , 
1  €  {1, 

Flow  sharing  problems  with  one  source  and  multiple  sinks  are  completely  analogous  to  the 
multiple- source  case:  merely  exchange  source  and  sinks  and  reverse  the  network.  For  the 
criteria  (ii)-(v),  one  can  even  allow  multiple  sources  and  multiple  sinks,  and  simultaneously 
optimize  one  criterion  for  the  sources  and  a  possibly  different  criterion  for  the  sinks.  This 
is  because  each  of  problems  (ii)-(v)  calls  for  a  maximum  flow,  and  a  multiple-source, 
multiple-sink  problem  can  be  decomposed  into  a  multiple-source  one-sink  problem  and 
a  one-source  multiple-sink  problem,  by  finding  a  minimum  cut  of  all  sources  from  all 
sinks,  contracting  all  vertices  on  the  sink  side  to  give  a  one-sink  problem,  and  separately 
contracting  all  vertices  on  the  source  side  to  give  a  one-source  problem.  This  observation 
is  due  to  Megiddo  [26]. 

Perfect  sharing  arises  in  a  network  transmission  problem  studied  by  Itai  and  Rodeh  [21]  and 
Gusfield  [17]  and  in  a  network  vulnerability  model  proposed  by  Cunningham  [5].  We  discuss 
these  models  below.  Brown  studied  maximin  sharing  [3],  Ichimori,  Ishii,  and  Nishida 
[20]  formulated  minimax  and  optimal  sharing,  and  Megiddo  [26,27]  studied  lexicographic 
sharing.  Motivation  for  these  problems  is  provided  by  the  following  kind  of  example,  which 
gives  rise  to  a  multiple  sink  problem.  During  a  famine,  relief  agencies  supplying  food  to  the 
stricken  areas  want  to  distribute  their  available  food  supplies  so  that  each  person  receives 
a  fair  share.  The  weight  associated  with  each  sink  (famine  area)  is  the  population  in 
that  area,  possibly  adjusted  for  differences  in  food  needs  between  adults  and  children,  and 
other  factors.  A  perfect  sharing  solution  gives  every  person  in  every  famine  area  the  same 
amount  of  food,  but  it  may  be  too  restrictive  since  it  need  not  allocate  all  the  available 
and  transportable  food  supply.  A  better  solution  will  be  provided  by  solving  one  of  the 
problems  (ii)-(v).  There  are  analogous  industrial  interpretations  of  this  model. 


We  shall  show  that  all  five  flow  sharing  problems  can  be  solved  in  0(nmlog(n2/m))  time 
using  the  algorithms  of  Sections  2  and  3.  The  lexicographic  sharing  problem  requires  com¬ 
puting  all  the  breakpoints  of  a  min-cut  capacity  function  by  the  algorithm  of  Section  3.3. 
The  other  four  problems  are  easier,  and  can  be  solved  by  the  algorithm  for  finding  the 
smallest  (or  largest)  breakpoint  given  in  Section  3.1.  Our  tool  for  solring  all  five  problems 
is  the  following  parametric  formulation:  for  each  Jt  £  5,  let  arc  (s,  st )  have  capacity  ji>fAf, 
where  A  is  a  real-valued  parameter.  Since  all  arc  capacities  are  nonnegative,  the  range  of 
interest  of  A  is  (0,oc).  There  are  at  most  k  breakpoints  of  the  min-cut  capacity  function 
k(A),  one  per  source  st. 

Perfect  sharing.  Find  the  smallest  breakpoint  A#  of  k(A).  Any  maximum  flow  for  Xt  solves 
the  perfect  sharing  problem.  This  was  observed  by  Gusfield  [17]  in  the  context  of  the 
network  transmission  scheduling  problem  described  below.  Another  application  will  arise 
in  Section  4.2. 

Scheduling  transmissions.  Itai  and  Rodeh  state  the  following  problem  of  scheduling  trans¬ 
missions  in  a  “circuit-switched”  communication  network  represented  by  a  directed  graph 
G  =  (V,  E)  with  fixed  positive  arc  capacities.  The  capacity  c(v,w)  is  the  effective  trans¬ 
mission  rate  of  the  communication  channel  (u,iu)  in  the  direction  from  v  to  w ,  say  in  bits 
per  second.  A  central  vertex  (sink)  t  £  V  receives  all  traffic  that  originates  at  a  subset 
of  vertices  5  C  V  -  {t}  called  emitters  (sources).  Each  emitter  s,  £  S  has  wt  >  0  bits  of 
information  that  it  wishes  to  send  to  t.  We  assume  that  G  has  paths  from  each  i,  £  5 
to  t.  The  communication  protocol  allows  the  sharing  of  arc  capacities  by  several  paths, 
but  it  requires  that  at  least  one  path  from  st  to  t  be  established  before  transmission  from 
can  begin.  Clearly,  if  transmissions  are  scheduled  from  each  emitter,  one  at  a  time, 
the  entire  task  can  be  completed  in  S»i€S  wi/e(XiiXi)  seconds,  where  (-Y,,A\)  is 

a  minimum  cut  separating  s ^  and  t.  But  since  arc  capacities  can  be  shared,  it  may  be 
possible  to  obtain  a  lower  value  for  T.  The  objective  is  to  minimize  the  time  T  within 
which  all  transmissions  can  be  completed. 

To  see  that  the  problem  is  a  perfect  sharing  multiple  source  problem,  let  A  =  1  /T,  in 
units  of  1/second,  and  assign  a  capacity  of  tt^A  bits  per  second  to  each  arc  (s,s,)  from  the 
supersource  s  to  an  emitter  £  S.  Once  A#  and  the  corresponding  maximum  flow  have 
been  computed  by  the  algorithm  described  above,  the  actual  transmission  schedule  can  be 
constructed  from  the  flow  in  O(m)  time  as  described  in  [21].  Itai  and  Rodeh  proposed  two 
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algorithms  for  this  problem,  with  running  times  of  0(kn2m)  and  0(k2nm log n).  These 
are  modifications  of  known  maximum  flow  algorithms.  In  comparison,  our  algorithm  runs 
in  O{nm\og(n2  /  m))  time. 

Maztmm  sharing.  Find  the  largest  breakpoint  Xt  of  /c(A).  Any  maximum  flow  for  Xt  solves 
the  maximin  sharing  problem. 

Mimmaz  sharing.  Find  the  smallest  breakpoint  A,  of  /c(A).  Find  a  maximum  flow  for 
A  .  Construct  a  residual  network  in  which  each  arc  (v,u’)  with  s  has  capacity 

c( v,w)  -  f(v,w),  each  arc  ( s ,  s{)  has  infinite  capacity,  and  each  arc  (st,  s)  has  zero  capacity. 
Find  a  maximum  flow  /'  in  the  residual  network.  The  flow  /  +  /'  is  a  minimax  flow  in  the 
original  network. 

Optimal  sharing.  Find  the  smallest  breakpoint  Xt  and  the  largest  breakpoint  A,  of  /c(A). 
Find  a  maximum  flow  /  for  Xt.  Construct  a  residual  network  in  which  each  arc  (r,  U’)  with 
s  #  {r.tr}  has  capacity  c(t,U’)  -  f(v,w)<  each  arc  (s,st)  has  capacity  u,(Aj  -  A4),  and 
each  arc  (^,4)  has  zero  capacity.  Find  a  maximum  flow  /'  in  the  residual  network.  The 
flow  f  ~  f'  is  an  optimal  flow  in  the  original  network. 

Lexicographic  sharing.  Find  all  the  breakpoints  of  k(A).  For  each  source  a,,  let  A,  be  the 
breakpoint  at  which  a,  moves  from  the  sink  side  to  the  source  side  of  a  minimum  cut.  For 
each  arc  (a,Sj)  define  its  capacity  to  be  tvtXi.  Find  a  maximum  flow  /  with  these  upper 
bounds  on  the  capacities  of  the  arcs  out  of  a.  Flow  /  is  a  lexicographic  flow  and  hence  an 
optimal  flow. 

The  correctness  of  the  first  four  algorithms  above  is  easy  to  verify.  We  shall  prove  the 
correctness  of  the  algorithm  for  the  lexicographic  sharing  problem.  Renumber  the  sources 
if  necessary  so  that  A}  <  A2  <  ■■•Xk,  and  let  G ^  denote  the  parametric  network  with 
A  =  oc. 

Theorem  4-1-  On  G ^  there  is  a  maximum  flow  /  such  that  /(a, a, )  =  tr.A,  for  all  i.  Such 
a  flow  is  a  lexicographic  flow. 

Proof.  Let  tj  ,t2,  •  •  • ,  i|_j  be  the  values  of  i  such  that  A,^  <  A^j.  Let  i0  =  0  and  i,  =  k. 
Then  At] ,  Atj,  •  •  • ,  Xti ,  are  the  distinct  breakpoints  in  increasing  order.  Let  {s  }  =  A’o  C 
A,  C  A’j  C  ■  •  •  C  A'j  be  the  sets  such  thai  (Xj,X})  for  1  <  j  <  l  is  the  minimum  cut 
with  the  smallest  sink  side  for  A  =  A,  .  Then  a,  .,,s,  G  A'-  -  A*  For 
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1  <  j  <  /,  the  cut  [Xj_1 , A’j.j)  is  a  minimum  cut  for  A  =  X{.  as  well,  specifically  the 
one  with  the  smallest  source  side.  Thus,  cAi  (A'j_1,A’j_1)  =  cA.  (Xj,X}).  It  follows  by 
induction  on  j  that  for  1  <  j  <  /, 


c> iixj>Xj)  =  Y,'t’iK+  Y  w*x'> 

i=]  i  =  » ,  +  1 


which  implies  that 


cooiXj  -  U),X})  =  YWiXii 


Equation  (4.2)  implies  that  any  flow  for  G ^  such  that  /(s,^)  =  wiXi  for  all  i  must  be  a 
maximum  flow  (choose  j  =  1  in  (4.2)).  It  must  also  be  a  lexicographic  flow,  since  for  all 
j,  any  flow  that  has  f(s,st)  >  il\ \i  for  1  <  i  <  t-_j  must  either  have  /(5,st)  =  in,  A,  for 
1  <  t  <  i:  or  have  some  t  £  {i j  -S-  1 ,  •  -  - ,  t  },  such  that  /(s,st  )  <  wt  A,. 

It  remains  to  show  that  Gx  admits  a  flow  /  with  /(s,Sj)  =  tn,Aj.  Consider  running  the 
min-cut  parametric  preflow  algorithm  presented  in  Section  2.5  on  the  parametric  network, 
for  the  successive  A  values  AM  ,  Atj ,  •  •  • ,  Atj .  Let  /j,/2,  •••,/(  be  the  successive  maximum 
preflows  generated  by  the  algorithm.  When  the  min-cut  preflow  algorithm  is  restarted  with 
a  new  value  A^  of  A,  the  flow  on  each  arc  (s,s,)  with  i  €  {*>-1  +!,••■  ,  k}  is  first  increased 
from  wt Xt)  s  to  u\ A^  .  All  of  this  new  flow  successfully  reaches  the  sink  t,  because  of 
equation  (4.1)  and  the  fact  that  (A';,  A';)  is  a  minimum  cut  for  A  =  A,  .  It  follows  by 
induction  on  j  that  is  a  flow  and  that  fj(s,st)  =  wiXl  for  1  <  t  <  t; .  In  particular,  ft 
is  the  desired  flow.  □ 

4.2.  Fractional  programming  applications 

Another  class  of  problems  that  can  be  solved  by  the  parametric  prefiow  algorithm  of 
Section  2.3  arises  from  various  discrete  and  network  optimization  problems  with  fractional 
objectives.  In  general,  the  fractional  programming  problem  is  defined  as 

A(r*)=  max  {  A(z)  =  f{x)/g{x)  }  ,  (4.3) 

r£S 

where  f{x)  and  g(x)  are  real-valued  functions  on  a  subset  5  of  J?n,  and  g(x)  >  0  for  all 
x  £  S .  Isbell  and  Marlow  |22  proposed  an  elegant  solution  method  for  the  important  case 
of  linear  /  and  g,  but  their  approach  has  been  extended  to  nonlinear  problems  (see  e.g. 


WWW 


Dinkelbach  [7]),  and  more  recently  to  several  classes  of  combinatorial  problems  (see  e.g. 
Picard  and  Queyranne  [32,  33],  Padberg  and  Wolsey  [30]  and  Cunningham  [5]). 

A  problem  that  is  intimately  related  to  (4.3)  is 


$ 
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r(x*,A)  =  max  {  z(x,A)  =  /(x)  -  A g(x)  }  ,  (4.4) 

*€S 

where  A  is  a  real- valued  constant.  These  two  problems  are  related  in  the  sense  that  x*  solves 
(4.3)  if  and  only  if  (x*,A*)  solves  (4.4)  for  A  =  A*  =  A(x*)  giving  the  value  z(x*,A*)  =  0. 
Isbell  and  Marlow’s  algorithm  generates  a  sequence  of  solutions  until  this  condition  is  met. 
We  state  their  algorithm  below  in  a  form  useful  for  our  purposes  (see  e.g.  Gondran  and 
Minoux  [14,  pp. 636-641]): 

Algorithm  FP: 

Step  0.  Select  some  x°  G  S.  Compute  A0  =  f(x°)/g(x°).  Set  k  =  0. 

Step  1.  Compute  xfc+1  solving  the  problem  (4.4):  z(xfc+1,Ak)  =  maxteS  z(x,  At)  = 

/(*)  - 

Step  2.  If  z(xfc+1,Ak)  =  0,  stop:  x*  =  x*.  Otherwise,  let  At+J  =  /(xt+I  )/</(x*+I ), 
replace  k  by  A:  -f  1  and  go  to  Step  1. 

Theorem  2.  Algorithm  FP  is  correct.  The  sequence  of  values  {Ak}  generated  by  the 
algorithm  is  decreasing. 

Proof.  For  any  particular  k,  z(x*+1,Afc)  is  nonnegative  in  Step  1,  since  z(x  k+\K)  > 
z(xfc,Ak)  =  0.  If  z(x*+1,Ak)  =  0,  the  algorithm  halts  with  xk'  which  solves  (4.4)  for 
A  =  Afc  and  hence  solves  (4.3).  The  algorithm  continues  only  if  z(xt+1,Ak)  >  0,  i.e. 
/(xfc+1)  —  A tg(xfc+1)  >  0,  which  implies  Afc+1  =  /(xfc+1  )/<?(xfc+I )  >  At.  □ 

If  maximization  is  replaced  by  minimization  in  problem  (4.3),  it  suffices  to  replace  max¬ 
imization  by  minimization  in  (4.4)  and  use  the  same  algorithm.  In  this  case  all  values 
of  z(x*'r,,Ak)  except  the  last  one  are  less  than  zero,  and  a  decreasing  sequence  {At}  is 
generated.  Another  important  observation  is  that  in  Step  1,  the  maximization  (4.4)  can 
be  taken  over  a  larger  set  S'  D  S,  provided  that  z(x,At)  <  0  for  all  x  G  S'  -  S  In  the 
minimization  problem,  the  corresponding  requirement  is  z(x,At)  >  0  for  all  x  G  S'  -  S. 
Several  of  our  applications  make  use  of  this  extension. 
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The  following  Lemma  can  be  used  to  bound  the  number  of  iterations  of  Algorithm  FP  in 
some  situations. 


Lemma  J.S.  </(x*+3)  <  9{xk)  for  k  >  1. 


Proof.  Consider  iterations  k  —  1  and  k  of  Algorithm  FP ,  and  assume  A(xfc)  <  A(z*).  In 
iteration  A  —  1  we  have  z(ar*,At_j)  >  0  and  Xk  =  f{xk)/ g(xk).  In  iteration  k  we  have: 

0  <  x(x*+1,A4)  =  /(x*+')-AlS(x‘+') 

=  /(x*+1)-  A4_!S(x*+1)  +  A4_l5(x*+1)  -  A4S(x*+1) 

<  /(**)  -  +  A,-,J(**+1)  -  A4j(*‘+') 

=  A4j(x*)  -  A4_,,(x*)  +  A4_l5(x**>)  -  A49(x*-') 

=  (s(x*)  -  s(x*+,))(At  -  A4_j), 

which  implies  that  g{xk)  >  p(ifc+1)  since  Xk  >  Xk_i.  □ 

The  efficiency  of  Algorithm  FP  depends  upon  the  number  of  times  problem  (4.4)  has  to 
be  solved,  and  on  the  time  spent  solving  it.  For  continuous  functions  f  and  g  defined  on 
a  nonempty  compact  set  S,  Schaible  [37]  has  shown  that  the  decreasing  sequence  {<7(1*)} 
for  k  >  1  approaches  g(x *)  linearly,  and  the  increasing  sequence  {At}  approaches  A* 
superlinearly.  Nevertheless,  (4.4)  may  be  as  hard  to  solve  as  the  original  fractional  problem 
unless  some  assumptions  are  made  about  /,  g  and  5.  Fortunately,  even  the  most  restrictive 
assumptions  find  relevant  applications  in  practice.  For  instance,  if  /  and  g  are  linear  and 
5  is  polyhedral  (the  case  in  [22]),  the  algorithm  consists  of  solving  a  finite  number  of 
linear  programs  (4.4)  whose  solution  is  implemented  by  cosi-parametric  programming  on 
intervals  [At,  At+ ,],  for  successive  k.  This  can  be  specialized  to  network  simplex  parametric 
programming  by  using  the  primitives  described  by  Grigoriadis  [15,.  If  /  is  a  negative 
6emidefinite  quadratic  form  and  g  is  linear,  the  sequence  of  concave  quadratic  programs 
defined  by  (4.4)  can  be  handled  by  the  parametric  algorithm  of  Grigoriadis  and  Ritter 
[16,.  If  /  and  g  are  negative  and  positive  definite  quadratic  forms  respectively,  Ritter’s 
parametric  quadratic  programming  method  [36]  can  be  used.  Approaches  for  more  general 
nonlinear  problems  are  analyzed  in  [7,  37j. 

If  S  is  nonempty  and  finite,  /  is  real-valued,  and  g  is  positive,  integer-valued,  and  bounded 
above  by  some  integer  p  >  0,  Lemma  4.3  implies  that  Algorithm  FP  will  terminate  in  p  + 1 
or  few’er  iterations.  This  observation  has  been  used  in  various  applications  where  p(x)  is  a 
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set  function,  for  which  usually  p  =  0(n).  Such  is  the  case  whether  (4.3)  is  a  maximization 
or  a  minimization  problem. 

We  shall  now  describe  a  number  of  applications  of  the  generic  Algorithm  FP.  In  each 
case  the  sequence  of  problems  (4.4)  that  arise  can  be  handled  by  our  parametric  preflow 
algorithm  or  its  min-cut  variant  described  in  Section  2.5. 

Strength  of  a  directed  network.  This  is  an  application  due  to  Cunningham  [5,  Section 
6].  Let  G  =  ( V,E )  be  a  given  directed  graph  with  n  vertices,  m  arcs,  nonnegative  arc 
weights  and  nonnegative  vertex  weights,  and  a  distinguished  vertex  s  £  V.  We  assume 
that  every  v  6  V  —  {s}  is  reachable  from  s  in  G.  The  arc  weight  c(v,tn)  represents 
the  cost  required  to  “destroy”  the  arc  (u,u>)  £  E.  The  node  weight  dv  is  the  “value” 
attributed  to  having  v  reachable  from  s.  Destroying  a  set  of  edges  A  C  E  (at  a  total 
cost  of  f(A)  =  ]C(t >  w)eAc(v'w)  )  ma.V  cause  some  subset  of  vertices  C  V  —  {s}  to 
become  unreachable  from  s,  resulting  in  a  loss  of  g{A)  =  52v€V  dv  in  total  value.  The 
ratio  f(A)/g(A)  is  the  cost  per  unit  reduction  in  value.  Cunningham  defines  the  strength 
of  the  network  to  be  the  minimum  of  this  ratio  taken  over  all  subsets  ACE  whose  removal 
reduces  the  value  of  the  network,  i.e.  such  that  g(A)  >  0.  This  is  a  problem  of  the  form 
(4.3): 

A(A')  =  min  {  A(4)  =  f(A)lg(A)  }  , 

ACE,9(A)>0 

which  leads  to  a  sequence  of  problems  (4.4)  that  Cunningham  calls  attack  problems: 
z(Afc+1,At)  =  max  {  *(A,Afc)  =  f(A)  -  Akg(A)  }  . 

ALE 

Each  such  problem  amounts  to  finding  a  minimum  cut  in  an  expanded  network  formed  by 
adding  to  G  a  sink  t  and  an  arc  (v,  t )  with  capacity  A kdv  for  each  v  £  V  — {s}.  If  we  solve  the 
strength  problem  using  Algorithm  FP  and  use  the  algorithm  of  Section  2.3  to  compute 
minimum  cuts  for  the  generated  parameter  values,  we  obtain  an  algorithm  running  in 
0(nmlog(n2/m))  time;  as  Cunningham  notes,  there  can  be  only  O(n)  iterations  of  Step 
1.  Alternatively,  we  can  make  use  of  his  observation  that  (4.5)  is  zero  if  and  only  if  there 
is  flow  in  the  expanded  network  such  that  f(v,t )  =  A kdv  for  each  v  £  V.  Equivalently, 
A(A*)  is  the  largest  value  of  A  for  which  the  minimum  cut  is  (!',  {<}).  That  is,  the  strength 
problem  is  a  perfect  sharing  problem,  and  it  can  be  solved  in  0(nmlog(n2/m))  time  as 
described  in  Section  4.1.  Either  method  improves  over  Cunningham’s  method,  which  solves 
0(n)  minimum  cut  problems  without  making  use  of  their  similarity. 
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Zero-one  fractional  programming.  An  important  subclass  of  (4.3)  is  the  problem  for  which 
f(x)  >  0  and  g(x)  >  0  are  given  polynomials  defined  for  all  x  in  5  =  {0,l}n  —  {0}n  as 
follows: 
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e  n  *• 

+  EX1* 

(4.5) 

pea  »€  p 

t=i 

s(x)  = 

£  n  ** 

+  x>. 

(4.6) 

Q€B  i€£? 

i=l 

The  sets  A  and  B  are  given  collections  of  nonempty  nonsingleton  subsets  of  {l,--*,n}, 
Op  >  0  for  each  P  €  A,  and  bq  <  0  for  each  Q  6  B.  Since  f(x)  >  0  and  g(x)  >  0  for  all 
x  €  5,  we  have  ai  >  0  and  b{  >  0  for  each  i  6  {1,  This  problem  was  studied  by 

Picard  and  Queyranne  [32].  For  ease  in  stating  time  bounds  we  assume  n  =  0(|A|  +  |B|). 

Algorithm  FP  leads  to  a  sequence  of  problems  of  the  form  (4.4)  for  increasing  values 
Xk  >  0  of  A.  Each  such  problem  is  an  instance  of  the  selection  or  provisioning  problem, 
characterized  by  Rhys  [35]  and  Balinski  [2]  as  a  minimum  cut  problem  on  a  bipartite  graph. 

The  entire  sequence  of  these  problems  can  be  handled  as  a  parametric  minimum  cut  prob¬ 
lem  of  the  kind  studied  in  Section  2.  We  give  two  different  formulations,  one  of  which 
works  for  the  special  case  of  B  =  0  (i.e.  g(x)  contains  no  nonlinear  terms)  and  the  other 
of  which  works  for  the  general  case.  All  the  subsequent  applications  we  consider  fall  into 
the  case  B  =  0. 

If  B  =  0,  we  define  a  bipartite  network  G  whose  vertex  set  contains  one  vertex  for  each 
set  P  €  A,  one  vertex  for  each  i  6  {1,**  • ,  n},  and  two  additional  vertices,  a  source  s  and 
a  sink  t.  There  is  an  arc  (s,r)  of  capacity  ap  for  each  vertex  v  corresponding  to  a  set 
P  €  A,  an  arc  (i,t)  of  capacity  Xbi  —  a,  for  every  i  6  {1,  •••,«},  and  an  arc  (v,i)  of  infinite 
capacity  for  every  vertex  corresponding  to  a  set  P  G  A  that  has  i  as  one  of  its  elements. 
Observe  that  the  capacities  of  all  arcs  into  t  are  nondecreasing  functions  of  A,  and  those 
of  all  oiher  arcs  are  constant.  The  parametric  preflow  algorithm  operates  on  GR  instead 
of  G.  Fo-  a  given  value  of  A,  a  minimum  cut  (X,  X)  in  G  corresponds  to  a  solution  x  to 
(4.4)  defined  by  xi  =  1  if  i  €  X,x{  —  0  if  t  €  X. 


In  the  general  case  ( B  ^  0),  it  is  convenient  to  assume  /(x)  >  0  for  some  x  (otherwise  the 
solution  to  (4.3)  is  A(x)  =  0,  attained  for  any  x)  and  that  algorithm  FP  starts  with  an 


*°  such  that  A0  >  0.  Then,  the  entire  sequence  {Xk}  is  strictly  positive.  To  solve  (4.4)  we 
rewrite  it  as  follows: 

z(x*,A)  =  A  max  (/(*)/ A  -  g(x))  . 

z£S 

We  define  the  network  G  to  have  a  vertex  set  consisting  of  one  vertex  for  each  set  P  6  A, 
one  vertex  for  each  set  Q  €  B,  one  vertex  for  each  t  6  {1,-  •  •  ,n},  and  a  source  s  and  a 
sink  t.  There  is  an  arc  (s,v)  of  capacity  aP/\  for  each  v  corresponding  to  a  set  P  £  A,  an 
arc  (s,v)  of  capacity  —  bq  for  each  v  corresponding  to  a  set  Q  £  B,  an  arc  (v,t)  of  infinite 
capacity  for  each  vertex  v  corresponding  to  a  set  P  £  A  or  Q  6  B  that  has  t  as  one  of  its 
elements,  and  an  arc  (*,<)  of  capacity  bi  —  ai/ A  for  each  t  £  {1, ••*,«}.  The  capacities  of 
arcs  out  of  s  are  nonincreasing  functions  of  A  and  those  of  arcs  into  t  are  nondecreasing 
functions  of  A;  the  parametric  preflow  algorithm  operates  on  GR.  Minimum  cuts  in  G 
correspond  to  solutions  exactly  as  described  above. 

Remark.  This  formulation  differs  from  that  in  [32]  because  of  the  division  by  A.  The 
formulation  of  [32]  gives  a  parametric  minimum  cut  problem  in  which  the  capacities  of 
arcs  out  of  the  source  and  of  arcs  into  the  sink  are  nondecreasing  functions  of  A,  to  which 
the  results  of  Section  2  do  not  apply. 

The  following  analysis  is  valid  for  both  of  the  above  network  formulations.  The  nesting 
property  of  minimum  cuts  (Lemma  2.4)  implies  that  the  number  of  iterations  of  Step  1  of 
Algorithm  FP  is  O(n),  a  fact  also  observed  by  Picard  and  Queyranne  [32].  To  state  time 
bounds,  let  us  denote  by  n'  and  m'  the  number  of  vertices  and  edges,  respectively,  in  G; 
n'  =  n  +  |A|  +  |I?|  +  2  and  m'  =  n  +  |A|  +  \B\  +  YIpza  +  X)c?€B  IQI-  Algorithm  FP,  in 
combination  with  the  parametric  preflow  algorithm  of  Section  2.3,  yields  a  time  bound  of 
0(n'm' log(n'2/m))  [or  0(nm'  \og(n2 /m'  +  2))],  improving  over  the  algorithms  of  Picard 
and  Queyranne  [32]  and  Gusfield,  Martel  and  Fernandez-Baca  [19]. 

Maximum-ratio  closure  problem.  This  problem  was  considered  by  Picard  and  Queyranne 
[33]  and  independently  by  Lawler  [24],  who  only  considered  acyclic  graphs  (see  the  next 
application).  The  problem  can  be  solved  by  a  straightforward  application  of  Algorithm 
FP.  Each  problem  in  the  sequence  of  problems  (4.4)  is  a  maximum- weight  closure  problem. 
The  maximum- weight  closure  problem  (Picard  [31])  is  the  generalization  to  nonbipbartite 
graphs  of  the  selection  or  provisioning  problem  [2,35]  mentioned  above. 


These  closure  problems  are  defined  formally  as  follows.  Let  G  =  {V,E)  be  a  directed 
graph  with  vertex  weights  av  of  arbitrary  sign.  A  subset  U  C  V  is  a  closure  in  G  if  for 
each  arc  (v,u>)  £  E  with  v  £  U  we  also  have  w  £  U.  A  closure  Um  C  V  is  of  maximum 
weight  if  the  sum  of  its  vertex  weights  is  maximum  among  all  closures  in  G.  To  compute  a 
maximum- weight  closure,  construct  the  graph  Gm  from  G  as  follows.  Add  a  source  s  and 
a  sink  t  to  G.  Create  an  arc  (s,t>)  of  capacity  av  and  an  arc  (v,t)  of  zero  capacity  for  each 
v  £  V.  Assign  infinite  capacity  to  all  arcs  in  E.  A  minimum  cut  (A”,Ar)  of  G *  gives  the 
desired  closure  Um  =  X  —  {a}. 

Now  let  av  >  0  and  bv  >  0  be  given  weights  on  the  vertices  of  G  =  (V',  E ).  The  maximum- 
ratio  closure  problem  is  to  find  a  closure  Um  that  maximizes  the  ratio  a{U)/b(U)  over  all 
nonempty  closures  17  C  V.  To  compute  a  maximum-ratio  closure,  Picard  and  Queyranne 
[33]  suggest  the  use  of  Algorithm  FP.  This  requires  the  solution  of  a  sequence  of  O(n) 
maximum-weight  closure  problems,  each  of  which  is  a  minimum  cut  problem.  Thus  an 
0(n2mlog(n2/m)-time  algorithm  results.  Lawler’s  algorithm  uses  binary  search  and  runs 
in  0(A:nmlog(n2/m))  time,  where  k  =  log(max{n, amax, &max}),  assuming  integer  weights. 

We  can  solve  the  entire  sequence  of  minimum  cut  problems  by  the  parametric  prefiow 
algorithm  of  Section  2.3  as  follows.  Modify  G *  so  that  for  each  vertex  v  £  V  there  is  an 
arc  (a,  v)  of  capacity  av  —  X bv  and  an  arc  (v,t)  of  capacity  zero.  All  other  arcs  have  infinite 
capacity.  We  start  with  U°  =  V  U  {a};  or,  equivalently,  with  a  sufficiently  small  value  of  A 
so  that  the  minimum  cut  is  ({a}  u  {<}).  Such  a  value  is  A0  =  min^  ajb i.  The  capacities 
of  arcs  out  of  the  source  are  nonincreasing  functions  of  the  parameter,  and  the  parameter 
values  are  given  on-line  in  increasing  order.  The  parametric  preflow  algorithm  operates  on 
and  runs  in  0(nmlog(n2/m))  time,  improving  the  bound  of  Picard  and  Queyranne 
by  a  factor  of  n  and  that  of  Lawler  by  a  factor  of  k. 

Remark.  Negative  arc  capacities  in  the  various  minimum  cut  problems  can  be  made  non¬ 
negative  using  the  transformation  suggested  in  Section  2.5.  In  the  minimum-ratio  closure 
problem,  it  suffices  to  assign  a  capacity  of  max{0,  av  -  Afe^}  to  each  arc  (s, v)  and  a  capacity 
of  n;<ix{0, Xbv  —  ae}  to  each  arc  (v,<). 

A  job  sequencing  application.  Lawler  [24]  applied  his  algorithm  to  a  problem  studied  by 
Sidney  [38]  and  others:  there  are  n  jobs  to  be  scheduled  for  processing  on  a  single  machine 
subject  to  a  partial  order  given  as  an  acyclic  graph  G  =  (V,  E),  where  V  is  the  set  of  jobs. 


Each  job  v  has  a  processing  time  av  and  a  “weight”  >  0  that  describes  its  importance  or 
some  measure  of  profit.  Let  the  completion  time  of  job  r  as  determined  by  a  given  feasible 
sequence  be  Cv.  It  is  required  to  find  a  sequence  that  minimizes  This  problem 

is  NP-complete  for  an  arbitrary  partial  order  even  when  all  av  =  1  or  all  bv  =  1  [24].  Sidney 
offered  the  following  decomposition  procedure.  First  find  a  maximum-ratio  closure  U1  such 
that  | £7,  |  is  minimum.  Remove  the  subgraph  induced  by  U1  from  G,  find  a  maximum-ratio 
closure  U2  of  the  reduced  graph,  and  repeat  this  process  until  the  entire  vertex  set  is 
partitioned.  Sidney  and  Lawler  call  closures  initial  sets  of  V.  Once  such  a  decomposition 
is  found,  an  optimal  schedule  can  be  computed  by  finding  an  optimal  schedule  for  each 
closure,  for  example  by  a  branch-and-bound  method,  and  then  concatenating  the  solutions. 
The  algorithm  described  above  can  be  used  to  find  each  closure.  The  overall  time  bound 
depends  upon  the  size  of  each  closure.  (Our  algorithm  will  give  closures  of  minimum 
cardinality,  since  the  algorithm  is  applied  to  the  graph  (G*)R;  see  Section  2.5.) 

Maximum  density  subgraph.  A  special  case  of  the  fractional  programming  problem  (4.3) 
is  that  of  finding  a  nonempty  subgraph  of  maximum  density  in  an  undirected  graph  G  = 
(V,  E )  with  n  vertices  and  m  edges.  The  density  of  a  subgraph  of  G  induced  by  a  subset  of 
vertices  V'  C  V'  is  the  number  of  its  edges  divided  by  the  number  of  its  vertices.  For  this 
application,  (4.5)  and  (4.6)  have  the  simpler  forms  f(x)  =  |xAx  and  g(x)  =  ex,  where  e 
is  the  vector  of  all  ones,  A  is  the  vertex- vertex  incidence  matrix  of  G,  and  xi  =  1  if  vertex 
i  €  T'  and  x^  =  0  otherwise.  Algorithm  FP  can  be  used  to  compute  a  maximum  density 
subgraph  of  G.  It  is  not  necessary  to  construct  a  bipartite  network  and  solve  minimum  cut 
problems  on  it.  We  can  merely  modify  G  by  specializing  the  construction  of  [34  .  Replace 
each  edge  of  G  by  two  oppositely  directed  arcs  of  unit  capacity,  add  a  source  s  and  a  sink 
t,  and  create  an  arc  (s,u)  of  capacity  6V  ~  X  and  an  arc  (v,f)  of  zero  capacity  for  each 
v  €  V-,  where  is  one- half  the  degree  of  vertex  v  in  G.  We  can  also  allow  weights  on  the 
edges  and  vertices.  The  resulting  algorithm  runs  in  G(nmlog(n2/m))  time.  This  bound 
is  better  than  that  of  Picard  and  Quevranne  [32]  by  a  factor  of  n,  and  better  than  that 
of  Goldberg  [12]  by  a  factor  of  log  n;  Goldberg’s  bound  is  valid  only  for  the  unweighted 
version  of  the  problem. 


4.3.  Parametric  zero-one  polynomial  functions 

We  consider  the  problem  of  computing  a  minimum  of  the  function 

z(A)  =  max  {  f(x )  -  A  (dr  -  b)  }  .  (4.7) 

*es 

where  5  =  {0,1}",  f(x)  is  a  polynomial  in  zero-one  variables  defined  by  (4.5),  and  di  >  0, 
i  €  {1,  •  •  •  ,n},  such  that  di  >  b  >  0. 

The  function  z(A)  differs  from  the  corresponding  function  (4.4)  that  arises  in  the  zero-one 
fractional  programming  application  of  Section  4.2  because  of  the  term  Xb  in  (4.7).  The 
function  z(A)  is  piecewise  linear  and  convex,  and  it  has  at  most  n  -  1  linear  segments 
and  n  —  2  breakpoints.  The  network  formulation  of  (4.7)  is  as  defined  for  the  zero-one 
fractional  programming  application.  The  breakpoints  of  z(A)  coincide  with  those  of  the 
min-cut  capacity  function  «(A)  for  this  network.  In  general,  no  minimum  of  z(A)  coincides 
w'ith  a  maximum  of  /c(A).  To  compute  a  minimum  of  z( A)  we  can  use  the  algorithm  of 
Section  3.2  for  finding  a  maximum  of  /c(A)  modified  to  use  the  graph  of  z(A)  instead  of  the 
graph  of  k(A)  to  guide  the  search.  We  have  z( 0)  =  o-p  +  at  >  0  (for  x  —  e). 

The  slope  of  the  leftmost  line  segment  of  z(A)  is  b  —  de  <  0,  and  the  slope  of  the  rightmost 
line  segment  is  6  >  0.  The  algorithm  consists  of  the  following  three  steps  and  finds  a 
minimum  of  z(A)  in  0(n'm' log(n,2/m'))  [or  0(nm'  log(n2/m'  -f  2))j  time.  A  cut  (A',  A') 
in  this  network  defines  a  solution  x  by  i,  =  1  if  vertex  t  €  A'  and  xt  =  0  otherwise. 

Step  0.  Start  with  initial  values  Aj  =  0,  =  e,  z(Aj)  =  /(a-1)  and  =  b  -  dxl . 

Choose  A 2  sufficiently  large  so  that  x2  =  0;  let  z(A2)  =  A26  and  f32  =  b. 

Step  1.  Compute  A3  =  (z(A2)  —  z(Aj  )/(0-l  —  /32)>  pass  A3  to  two  invocations,  I  and  D, 
of  the  parametric  preflow  algorithm,  and  run  them  concurrently.  If  invocation  I  finds 
a  minimum  cut  (Xf,X2)  first,  suspend  invocation  D  and  go  to  Step  2  (the  other  case 
is  symmetric.)  Compute  (3l  =  b  —  dx3 ,  the  slope  of  the  line  segment  of  z(A)  derived 
from  this  cut. 

Step  t.  If  03  =  0,  stop:  A*  =  A3.  Otherwise,  if  03  >  0,  replace  A2  by  A3,  back  up 
invocation  D  to  its  state  before  it  began  processing  A3,  and  go  to  Step  1.  Otherwise, 
go  to  Step  3. 

Step  3  (03  <  0).  Finish  the  invocation  D  for  A3.  This  produces  a  minimum  cut 
(XP^P),  not  necessarily  the  same  as  (Xj,X3).  If  03  >  0,  stop:  A*  =  A3.  Otherwise, 


replace  Aj  by  As,  back  up  invocation  I  to  its  state  before  processing  As,  and  go  to 
Step  1. 

We  now  describe  an  application  of  this  algorithm. 

Knapsack-constrained  provisioning  problems.  We  consider  the  following  provisioning  prob¬ 
lem  with  a  knapsack  constraint  that  limits  the  weight  of  the  selected  items: 

z{xm)  =  *  max  ^  {  /( x)  |  dx  <  b  },  (4.8) 

where  /( i)  is  given  by  (4.5),  d  is  a  positive  n-vector  of  item  weights  and  b  is  a  scalar,  the 
knapsack  size,  such  that  J2iev  d,  >  b  >  0,  V  =  {!,-  •  •  ,n}.  Thus,  in  addition  to  the  benefit 
ai  obtained  for  including  an  individual  item  t  £  l’  in  the  knapsack,  the  model  allows  the 
possibility  of  an  additional  reward  of  ap  >  0  for  including  all  of  the  items  that  comprise 
a  given  subset  P  €  A.  The  (linear)  knapsack  problem  is  a  special  case  of  (4.8)  in  which  all 
subsets  P  €  A  are  singletons. 

This  NP-complete  problem  was  suggested  by  Lawler  [23j.  Because  of  its  many  practical 
applications  there  is  interest  in  the  fast  computation  of  bounds  on  z(r’).  To  this  end  we 
consider  the  Lagrangean  function  for  (4.8): 

L(x,  A)  =  f(x)  —  A  (dx  -  b ),  for  A  >  0, 

which  has  a  finite  infimum  over  x  6  {0,  l}n.  For  each  A  >  0,  we  define  the  dual  function: 

$(A)  =  max  L(x,  A). 

*€{0,1}" 

$(A)  is  a  piecewise  linear  convex  function  of  A,  having  at  most  n  —  1  line  segments  and 
n  —  2  breakpoints.  We  wish  to  solve  the  following  Lagrange  an  dual  problem: 

$(A*)  =  min  $(A). 

A>0 

This  value  is  an  upper  bound  on  z(z*)  and  can  be  used  to  construct  heuristics  and  search 
procedures  for  computing  an  approximate  or  exact  solution  to  (4.8).  It  can  be  evaluated 
by  the  above  algorithm  in  0(n'm'  log(n'2/m'))  (or  0(nm'  log(n2/m'  +  2))]  time. 

A  special  case  of  considerable  practical  importance  is  the  quadratic  knapsack  problem,  for 
which  /(x)  =  xAx  where  A  =  [a^J  is  a  nonnegative  real  symmetric  matrix  having  no 


null  rows.  For  this  case,  Gallo,  Hammer  and  Simeone  [11]  proposed  an  0(n2  logn)-time 
algorithm  for  creating  a  class  of  “upper  planes”  bounding  z(z).  Chaillou,  Hansen  and 
Mahieu  [4]  showed  that  its  Lagrangean  dual  can  be  solved  as  a  sequence  of  0(n)  minimum 
cut  problems  in  0(n2mlog(n2/m))  time. 

The  problem  of  evaluating  $(A)  for  a  fixed  A  can  be  formulated  as  a  minimum  cut  problem 
using  a  graph  construction  similar  to  that  described  earlier  for  the  maximum  density 
subgraph  problem,  thereby  avoiding  the  use  of  a  bipartite  graph.  Let  G  =  (V,  E)  be  a 
directed  graph  with  vertex  set  V  =  {l, •••,»},  arc  set  E  =  {(r,tn)  :  aeu,  >  0,v,u>  £  V}, 
and  arc  weights  a(v,w)  =  avw.  We  add  to  G  a  source  s,  a  sink  t,  and  an  arc  (i,u)  of 
capacity  av  —  A dv  and  an  arc  (v,t)  of  zero  capacity  for  each  v  £  V.  Using  this  network 
formulation,  the  above  algorithm  computes  the  Lagrangean  relaxation  of  a  quadratric 
knapsack  problem  in  0(nmlog(n2/m))  time. 

4.4.  Miscellaneous  applications 

Our  last  two  applications  both  use  the  algorithm  developed  in  Section  3.3  for  computing  the 
min-cut  capacity  function  *(A)  of  a  parametric  minimum  cut  problem.  The  first  application 
is  to  a  problem  of  computing  critical  load  factors  for  modules  of  a  distributed  program  in 
a  two-processor  distributed  system  [42] .  The  second  application  is  to  a  problem  of  record 
segmentation  between  primary  and  secondary  memory  in  large  shared  databases  [9]. 

Critical  load  factors  in  two -process or  distributed  systems.  Stone  [42]  modeled  this  problem 
by  a  graph  G  =  (V,  E)  in  which  V  =  is  the  set  of  program  modules  and  E  is 

the  set  of  pairs  of  modules  that  need  to  communicate  with  each  other.  The  capacity  of  an 
arc  (v,tr)  £  E  specifies  the  communication  cost  between  modules  v  and  w  (it  is  infinity 
if  the  modules  must  be  coresident).  The  two  processors,  say  A  and  B,  are  represented  by 
the  source  s  and  the  sink  t ,  respectively,  that  are  appended  to  the  network.  There  is  an 
arc  (a,v)  of  capacity  A bv  >  0  where  6„  is  the  given  cost  of  executing  program  module  v  on 
processor  B.  There  is  an  arc  (v,t)  of  capacity  (1  -  A )av  >  0  where  av  is  the  given  cost  of 
executing  program  module  v  on  processor  A. 

The  parameter  A  £  [0,1]  is  the  fraction  of  the  time  processor  A  delivers  useful  cycles, 
commonly  known  as  the  load  factor.  For  a  fixed  value  of  A,  a  minimum  cut  (A',  A')  in  this 
network  gives  an  optimum  assignment  of  modules  to  processors.  For  A  =  0,  a  minimum 
cut  (X,X)  with  jXj  minimum  has  X  =  {s},  i.e.  all  modules  are  assigned  to  B.  For  A  =  1 


a  minimum  cut  (X,X)  with  |Jf|  maximum  has  X  =  {<},  i.e.  all  modules  are  assigned 
to  A.  The  objective  is  to  find  the  best  assignment  of  program  modules  to  processors  for 
various  values  of  A,  or  to  generate  these  assignments  for  each  breakpoint  of  the  min-cut 
capacity  function  k(A).  Lemma  2.4  implies  that,  at  each  breakpoint,  one  or  more  modules 
shift  from  one  side  of  the  cut  to  the  other.  By  listing,  for  each  module,  the  breakpoint 
at  which  it  shifts  from  one  side  of  the  minimum  cut  to  the  other,  one  can  determine 
what  Stone  calls  the  critical  load  factor  for  each  module.  The  operating  system  can  then 
use  this  list  of  critical  load  factors  to  do  dynamic  assignment  of  modules  to  processors. 
The  algorithm  of  Section  3.3  will  compute  the  critical  load  factors  of  the  modules  in 
0(nm'log((n  +  2 )2/m'))  time,  where  m!  =  m  +  2n.  Stone  does  not  actually  propose  an 
algorithm  for  this  computation. 

Record  segmentation  in  large  shared  databases.  Eisner  and  Severance  [9  have  stated  a 
model  for  segmenting  records  in  a  large  shared  database  between  primary  and  secondary 
memory.  Such  a  database  consists  of  a  set  of  data  items  S  =  {1,  •  •  • ,  N]  and  serves  a  set 
of  users  T  =  {1,  •••,«}.  Each  user  u>  G  T  retrieves  a  nonempty  subset  Su.  C  S  of  data 
items  and  receives  a  “value”  (satisfaction)  of  hu.  >  0  whenever  all  of  the  items  in  Sw  reside 
in  primary  memory.  The  cost  of  transporting  and  storing  a  data  item  r  €  S  in  primary 
memory  is  A av  >  0,  where  a„  >  0.  The  scalar  A  >  0  is  a  conversion  factor  such  that  A 
units  of  transportation  and  storage  costs  equals  one  unit  of  user  value.  The  objective  is  to 
find  a  segmentation  that  minimizes  the  total  cost  minus  user  satisfaction. 


For  a  fixed  value  of  A  the  problem  can  be  formulated  as  a  selection  or  provisioning  problem 
[2,35  as  follows.  Construct  a  bipartite  graph  having  the  data  items  5  as  its  source  part 
and  the  users  T  as  its  sink  part.  Construct  an  arc  (v,w),  v  €  5,  iv  €  T  of  infinite  capacity 
if  data  item  v  belongs  to  the  set  of  data  items  Sw  retrieved  by  user  w.  Create  a  supersource 
s  and  a  supersink  t ,  and  append  an  arc  (s,v)  of  capacity  Xav  for  each  r  £  5  and  an  arc 
(w,t)  of  capacity  for  each  w  £  T.  A  min-cut  (X, X)  separating  s  and  i  in  this  network 
necessarily  partitions  S  and  T  into  {SX,SX)  and  (TX,TX),  respectively.  It  is  easy  to  see 
that 


c{X,X)  =  min  -v 

SxuTx,SxuTx 

v€Sx 


{ *  y.  *„+  y.  >• 

w£Tx 


The  value  of  A  plays  an  important  role  in  this  linear  performance  measure,  and  it  depends 
upon  the  system  load.  In  practice  it  is  necessary  to  create  a  list  of  primary  storage  assign¬ 
ments  for  all  critical  values  of  A.  The  database  inquiry  program  can  then  select  and  imple- 


xnent  the  best  assignment  at  appropriate  times.  This  table  consists  of  all  the  breakpoints  of 
the  min-cut  capacity  function  k(A)  and,  for  each  data  item  and  user,  the  parameter  value 
at  which  it  moves  from  one  side  to  the  other  of  a  minimum  cut.  This  information  can  be 
computed  by  the  breakpoint  algorithm  of  Section  3.3  in  0((n  +  lV)mlog((n  +  N)2 /m))  [or 
0(min{n,iST}mlog((min{n,  JV})2/m  +  2))J  time.  The  algorithm  proposed  by  Eisner  and 
Severance  for  solving  the  parametric  problem  requires  the  solution  of  0(min{n,  N})  mini¬ 
mum  cut  problems.  Our  algorithm  improves  their  method  by  a  factor  of  min{n,  N}.  They 
also  consider  a  nonlinear  performance  measure,  for  which  an  algorithm  like  that  in  Section 
4.3  can  be  used  to  derive  bounds  on  an  optimum  solution.  This  bounding  method  gives 
an  approximate  solution,  and  the  method  can  be  used  in  a  branch-and-bound  algorithm 
to  give  an  exact  solution. 

5.  Remarks 

We  have  shown  how  to  extend  the  maximum  flow  algorithm  of  Goldberg  and  Tarjan  to 
solve  a  sequence  of  0(n)  related  maximum  flow  problems  at  a  cost  of  only  a  constant 
factor  over  the  time  to  solve  one  problem.  The  problems  must  be  instances  of  the  same 
parametric  maximum  flow  problem  and  the  corresponding  parameter  values  must  either 
consistently  increase  or  consistently  decrease.  We  have  further  shown  how  to  extend  the 
algorithm  to  generate  the  entire  min-cut  capacity  function  of  such  a  parametric  problem, 
assuming  that  the  arc  capacities  are  linear  functions  of  the  parameter. 

We  have  applied  our  algorithms  to  solve  a  variety  of  combinatorial  optimization  problems, 
deriving  improved  time  bounds  for  each  of  the  problems  considered.  Our  list  of  applications 
is  meant  to  be  illustrative,  not  exhaustive.  WTe  expect  that  more  applications  will  be 
discovered.  Although  we  have  only  considered  a  special  form  of  the  parametric  maximum 
flow  problem,  most  of  the  parametric  maximum  flow’  problems  we  have  encountered  in  the 
literature  can  be  put  into  this  special  form. 

We  have  discussed  only  sequential  algorithms  in  this  paper,  but  our  ideas  extend  to  the 
realm  of  parallel  algorithms.  Specifically,  the  preflow  algorithm  has  a  parallel  version 
that  runs  in  0(n2  log  n)  time  using  n  processors  on  a  parallel  random-access  machine. 
This  version  extends  to  the  parametric  preflow  algorithm  in  exactly  the  same  wrav  as  the 
sequential  algorithm.  Thus  we  obtain  0(n2  logn)-time,  n-processor  parallel  algorithms  for 


the  problems  considered  in  Sections  2  and  3  and  for  each  of  the  applications  in  Section  4, 
where  n  is  the  number  of  vertices  in  the  network. 


There  are  a  number  of  remaining  open  problems.  One  is  to  find  additional  applications. 
Possible  applications  include  computing  the  arboricity  of  a  graph  [29,32]  and  computing 
properties  of  activity  selection  games  [44].  Gusfield  [18  has  recently  found  a  new  ap¬ 
plication,  to  a  problem  considered  by  Cunningham  [5],  of  solving  the  sequence  of  attack 
problems  involved  in  the  computation  of  the  strength  of  an  undirected  graph.  (This  prob¬ 
lem  is  related  to  the  strength  problem  considered  in  Section  4.2  but  is  harder.) 

Another  area  for  research  is  investigating  whether  an  arbitrary  maximum  flow  algorithm 
can  be  extended  to  the  parametric  problem  at  a  cost  of  only  a  constant  factor  in  running 
time.  One  algorithm  that  we  have  unsuccessfully  tried  to  extend  in  this  way  is  that  of 
Ahuja  and  Orlin  [l] .  Working  in  this  direction,  Martel  [25]  has  recently  discovered  how  to 
modify  an  algorithm  based  on  the  approach  of  Dinic  f6’  so  that  it  solves  the  parametric 
problem  with  only  a  constant  factor  increase. 
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